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

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

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

revision 1.8 by pam-fi, Wed Oct 25 16:18:41 2006 UTC revision 1.31 by pam-fi, Fri Aug 31 14:56:52 2007 UTC
# Line 12  Line 12 
12        subroutine track_finding(iflag)        subroutine track_finding(iflag)
13    
14        include 'commontracker.f'        include 'commontracker.f'
15          include 'level1.f'
16        include 'common_momanhough.f'        include 'common_momanhough.f'
17        include 'common_mech.f'        include 'common_mech.f'
18        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
19        include 'common_mini_2.f'        include 'common_mini_2.f'
20        include 'calib.f'        include 'calib.f'
       include 'level1.f'  
21        include 'level2.f'        include 'level2.f'
22    
       include 'momanhough_init.f'  
         
 c      logical DEBUG  
 c      common/dbg/DEBUG  
23    
24    
25    c      print*,'======================================================'
26    c$$$      do ic=1,NCLSTR1
27    c$$$         if(.false.
28    c$$$     $        .or.nsatstrips(ic).gt.0
29    c$$$c     $        .or.nbadstrips(0,ic).gt.0
30    c$$$c     $        .or.nbadstrips(4,ic).gt.0
31    c$$$c     $        .or.nbadstrips(3,ic).gt.0
32    c$$$     $        .or..false.)then
33    c$$$            print*,'--- cl-',ic,' ------------------------'
34    c$$$            istart = INDSTART(IC)
35    c$$$            istop  = TOTCLLENGTH
36    c$$$            if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
37    c$$$            print*,'ADC   ',(CLADC(i),i=istart,istop)
38    c$$$            print*,'s/n   ',(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
39    c$$$            print*,'sgnl  ',(CLSIGNAL(i),i=istart,istop)
40    c$$$            print*,'strip ',(i-INDMAX(ic),i=istart,istop)
41    c$$$            print*,'view ',VIEW(ic)
42    c$$$            print*,'maxs ',MAXS(ic)
43    c$$$            print*,'COG4 ',cog(4,ic)
44    c$$$            ff = fbad_cog(4,ic)
45    c$$$            print*,'fbad ',ff
46    c$$$            print*,(CLBAD(i),i=istart,istop)
47    c$$$            bb=nbadstrips(0,ic)
48    c$$$            print*,'#BAD (tot)',bb
49    c$$$            bb=nbadstrips(4,ic)
50    c$$$            print*,'#BAD (4)',bb
51    c$$$            bb=nbadstrips(3,ic)
52    c$$$            print*,'#BAD (3)',bb
53    c$$$            ss=nsatstrips(ic)
54    c$$$            print*,'#saturated ',ss
55    c$$$         endif
56    c$$$      enddo
57          
58  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
59  *     STEP 1  *     STEP 1
60  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
# Line 44  c      common/dbg/DEBUG Line 74  c      common/dbg/DEBUG
74  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
75  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
76    
77  c      iflag=0  
78        call cl_to_couples(iflag)        call cl_to_couples(iflag)
79        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
80           goto 880               !fill ntp and go to next event                       goto 880               !go to next event
81        endif        endif
82                
 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  
 *     selezione di tracce pulite per diagnostica  
 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  
 c$$$         if(DEBUG)then  
 c$$$            do ip=1,nplanes  
 c$$$               if(ncp_plane(ip).ne.1)good2=.false.  
 c$$$            enddo  
 c$$$c            if(good2.eq.0)goto 100!next event  
 c$$$c            if(good2.eq.0)goto 880!fill ntp and go to next event  
 c$$$         endif  
           
   
   
83  *-----------------------------------------------------  *-----------------------------------------------------
84  *-----------------------------------------------------  *-----------------------------------------------------
85  *     HOUGH TRASFORM  *     HOUGH TRASFORM
# Line 91  c$$$         endif Line 108  c$$$         endif
108  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
109  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
110    
111  c      iflag=0  
112        call cp_to_doubtrip(iflag)        call cp_to_doubtrip(iflag)
113        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
114           goto 880               !fill ntp and go to next event                       goto 880               !go to next event            
115        endif        endif
116                
117                
# Line 123  c      iflag=0 Line 140  c      iflag=0
140  *     $     ,ptcloud_xz,tr_cloud,cpcloud_xz            *     $     ,ptcloud_xz,tr_cloud,cpcloud_xz          
141  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
142  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
143    *     count number of hit planes
144          planehit=0                
145          do np=1,nplanes          
146            if(ncp_plane(np).ne.0)then  
147              planehit=planehit+1  
148            endif                  
149          enddo                    
150          if(planehit.lt.3) goto 880 ! exit              
151          
152          nptxz_min=x_min_start              
153          nplxz_min=x_min_start              
154                
155          nptyz_min=y_min_start              
156          nplyz_min=y_min_start              
157    
158          cutdistyz=cutystart      
159          cutdistxz=cutxstart      
160    
161  c      iflag=0   878  continue
162        call doub_to_YZcloud(iflag)        call doub_to_YZcloud(iflag)
163        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
164           goto 880               !fill ntp and go to next event                       goto 880               !fill ntp and go to next event            
165        endif        endif
166  c      iflag=0        if(nclouds_yz.eq.0.and.cutdistyz.lt.maxcuty)then
167            if(cutdistyz.lt.maxcuty/2)then
168              cutdistyz=cutdistyz+cutystep
169            else
170              cutdistyz=cutdistyz+(3*cutystep)
171            endif
172            goto 878                
173          endif                    
174    
175          if(planehit.eq.3) goto 881          
176          
177     879  continue  
178        call trip_to_XZcloud(iflag)        call trip_to_XZcloud(iflag)
179        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
180           goto 880               !fill ntp and go to next event                       goto 880               !fill ntp and go to next event            
181        endif        endif
182                                    
183          if(nclouds_xz.eq.0.and.cutdistxz.lt.maxcutx)then
184            cutdistxz=cutdistxz+cutxstep
185            goto 879                
186          endif                    
187    
188        
189     881  continue  
190    *     if there is at least three planes on the Y view decreases cuts on X view
191          if(nclouds_xz.eq.0.and.nclouds_yz.gt.0.and.
192         $     nplxz_min.ne.y_min_start)then
193            nptxz_min=x_min_step    
194            nplxz_min=x_min_start-x_min_step              
195            goto 879                
196          endif                    
197            
198   880  return   880  return
199        end        end
200    
# Line 144  c      iflag=0 Line 204  c      iflag=0
204        subroutine track_fitting(iflag)        subroutine track_fitting(iflag)
205    
206        include 'commontracker.f'        include 'commontracker.f'
207          include 'level1.f'
208        include 'common_momanhough.f'        include 'common_momanhough.f'
209        include 'common_mech.f'        include 'common_mech.f'
210        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
211        include 'common_mini_2.f'        include 'common_mini_2.f'
212        include 'calib.f'        include 'calib.f'
       include 'level1.f'  
213        include 'level2.f'        include 'level2.f'
214    
215        include 'momanhough_init.f'  c      include 'momanhough_init.f'
216                
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
217        logical FIMAGE            !        logical FIMAGE            !
218          real trackimage(NTRACKSMAX)
219          real*8 AL_GUESS(5)
220    
221  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
222  *     STEP 4   (ITERATED until any other physical track isn't found)  *     STEP 4   (ITERATED until any other physical track isn't found)
# Line 198  c         iflag=0 Line 257  c         iflag=0
257           ibest=0                !best track among candidates           ibest=0                !best track among candidates
258           iimage=0               !track image           iimage=0               !track image
259  *     ------------- select the best track -------------  *     ------------- select the best track -------------
260           rchi2best=1000000000.  c$$$         rchi2best=1000000000.
261    c$$$         do i=1,ntracks
262    c$$$            if(RCHI2_STORE(i).lt.rchi2best.and.
263    c$$$     $         RCHI2_STORE(i).gt.0)then
264    c$$$               ibest=i
265    c$$$               rchi2best=RCHI2_STORE(i)
266    c$$$            endif
267    c$$$         enddo
268    c$$$         if(ibest.eq.0)goto 880 !>> no good candidates
269    
270    *     -------------------------------------------------------
271    *     order track-candidates according to:
272    *     1st) decreasing n.points
273    *     2nd) increasing chi**2
274    *     -------------------------------------------------------
275             rchi2best=1000000000.
276             ndofbest=0            
277           do i=1,ntracks           do i=1,ntracks
278              if(RCHI2_STORE(i).lt.rchi2best.and.             ndof=0              
279       $         RCHI2_STORE(i).gt.0)then             do ii=1,nplanes    
280                 ndof=ndof        
281         $            +int(xgood_store(ii,i))
282         $            +int(ygood_store(ii,i))
283               enddo              
284               if(ndof.gt.ndofbest)then
285                 ibest=i
286                 rchi2best=RCHI2_STORE(i)
287                 ndofbest=ndof    
288               elseif(ndof.eq.ndofbest)then
289                 if(RCHI2_STORE(i).lt.rchi2best.and.
290         $            RCHI2_STORE(i).gt.0)then
291                 ibest=i                 ibest=i
292                 rchi2best=RCHI2_STORE(i)                 rchi2best=RCHI2_STORE(i)
293              endif                 ndofbest=ndof  
294           enddo               endif            
295               endif
296             enddo
297    
298    c$$$         rchi2best=1000000000.
299    c$$$         ndofbest=0             !(1)
300    c$$$         do i=1,ntracks
301    c$$$           if(RCHI2_STORE(i).lt.rchi2best.and.
302    c$$$     $          RCHI2_STORE(i).gt.0)then
303    c$$$             ndof=0             !(1)
304    c$$$             do ii=1,nplanes    !(1)
305    c$$$               ndof=ndof        !(1)
306    c$$$     $              +int(xgood_store(ii,i)) !(1)
307    c$$$     $              +int(ygood_store(ii,i)) !(1)
308    c$$$             enddo              !(1)
309    c$$$             if(ndof.ge.ndofbest)then !(1)
310    c$$$               ibest=i
311    c$$$               rchi2best=RCHI2_STORE(i)
312    c$$$               ndofbest=ndof    !(1)
313    c$$$             endif              !(1)
314    c$$$           endif
315    c$$$         enddo
316    
317           if(ibest.eq.0)goto 880 !>> no good candidates           if(ibest.eq.0)goto 880 !>> no good candidates
318  *-------------------------------------------------------------------------------      *-------------------------------------------------------------------------------    
319  *     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 332  c         iflag=0
332              iimage=0              iimage=0
333           endif           endif
334           if(icand.eq.0)then           if(icand.eq.0)then
335              print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand              if(VERBOSE.EQ.1)then
336       $           ,ibest,iimage                 print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand
337         $              ,ibest,iimage
338                endif
339              return              return
340           endif           endif
341    
# Line 236  c         iflag=0 Line 346  c         iflag=0
346  *     **********************************************************  *     **********************************************************
347  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
348  *     **********************************************************  *     **********************************************************
349             call guess()
350             do i=1,5
351                AL_GUESS(i)=AL(i)
352             enddo
353    c         print*,'## guess: ',al
354    
355           do i=1,5           do i=1,5
356              AL(i)=dble(AL_STORE(i,icand))                          AL(i)=dble(AL_STORE(i,icand))            
357           enddo           enddo
358            
359           IDCAND = icand         !fitted track-candidate           IDCAND = icand         !fitted track-candidate
360           ifail=0                !error flag in chi2 computation           ifail=0                !error flag in chi2 computation
361           jstep=0                !# minimization steps           jstep=0                !# minimization steps
362    
363           iprint=0           iprint=0
364           if(DEBUG)iprint=1  c         if(DEBUG.EQ.1)iprint=1
365             if(VERBOSE.EQ.1)iprint=1
366             if(DEBUG.EQ.1)iprint=2
367           call mini2(jstep,ifail,iprint)           call mini2(jstep,ifail,iprint)
368           if(ifail.ne.0) then           if(ifail.ne.0) then
369              if(DEBUG)then              if(VERBOSE.EQ.1)then
370                 print *,                 print *,
371       $              '*** MINIMIZATION FAILURE *** (mini2) '       $              '*** MINIMIZATION FAILURE *** (after refinement) '
372       $              ,iev       $              ,iev
373    
374    c$$$               print*,'guess:   ',(al_guess(i),i=1,5)
375    c$$$               print*,'previous: ',(al_store(i,icand),i=1,5)
376    c$$$               print*,'result:   ',(al(i),i=1,5)
377    c$$$               print*,'xgood ',xgood
378    c$$$               print*,'ygood ',ygood
379    c$$$               print*,'----------------------------------------------'
380              endif              endif
381              chi2=-chi2  c            chi2=-chi2
382           endif           endif
383                    
384           if(DEBUG)then           if(DEBUG.EQ.1)then
385              print*,'----------------------------- improved track coord'              print*,'----------------------------- improved track coord'
386  22222       format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)  22222       format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
387              do ip=1,6              do ip=1,6
# Line 265  c         iflag=0 Line 392  c         iflag=0
392           endif           endif
393    
394  c         rchi2=chi2/dble(ndof)  c         rchi2=chi2/dble(ndof)
395           if(DEBUG)then           if(DEBUG.EQ.1)then
396              print*,' '              print*,' '
397              print*,'****** SELECTED TRACK *************'              print*,'****** SELECTED TRACK *************'
398              print*,'#         R. chi2        RIG'              print*,'#         R. chi2        RIG'
# Line 281  c         rchi2=chi2/dble(ndof) Line 408  c         rchi2=chi2/dble(ndof)
408  *     ------------- search if the track has an IMAGE -------------  *     ------------- search if the track has an IMAGE -------------
409  *     ------------- (also this is stored )           -------------  *     ------------- (also this is stored )           -------------
410           if(FIMAGE)goto 122     !>>> jump! (this is already an image)           if(FIMAGE)goto 122     !>>> jump! (this is already an image)
411  *     now search for track-image, by comparing couples IDs  
412    *     -----------------------------------------------------
413    *     first check if the track is ambiguous
414    *     -----------------------------------------------------
415    *     (modified on august 2007 by ElenaV)
416             is1=0
417             do ip=1,NPLANES
418                if(SENSOR_STORE(ip,icand).ne.0)then
419                   is1=SENSOR_STORE(ip,icand)
420                   if(ip.eq.6)is1=3-is1 !last plane inverted
421                endif
422             enddo
423             if(is1.eq.0)then
424                if(WARNING.EQ.1)print*,'** WARNING ** sensor=0'
425                goto 122            !jump
426             endif
427    c         print*,'is1 ',is1
428             do ip=1,NPLANES
429                is2 = SENSOR_STORE(ip,icand) !sensor
430    c            print*,'is2 ',is2,' ip ',ip
431                if(ip.eq.6.and.is2.ne.0)is2=3-is2 !last plane inverted
432                if(
433         $           (is1.ne.is2.and.is2.ne.0)
434         $           )goto 122      !jump (this track cannot have an image)
435             enddo
436             if(DEBUG.eq.1)print*,' >>> ambiguous track! '
437    *     now search for track-image among track candidates
438    c$$$         do i=1,ntracks
439    c$$$            iimage=i
440    c$$$            do ip=1,nplanes
441    c$$$               if(     CP_STORE(nplanes-ip+1,icand).ne.
442    c$$$     $              -1*CP_STORE(nplanes-ip+1,i).and.
443    c$$$     $              CP_STORE(nplanes-ip+1,i).ne.0.and.
444    c$$$     $              CP_STORE(nplanes-ip+1,icand).ne.0 )iimage=0
445    c$$$               print*,' track ',i,' CP ',CP_STORE(nplanes-ip+1,i)
446    c$$$     $              ,CP_STORE(nplanes-ip+1,icand),iimage
447    c$$$            enddo
448    c$$$            if(  iimage.ne.0.and.
449    c$$$c     $           RCHI2_STORE(i).le.CHI2MAX.and.
450    c$$$c     $           RCHI2_STORE(i).gt.0.and.
451    c$$$     $           .true.)then
452    c$$$               if(DEBUG.EQ.1)print*,'Track candidate ',iimage
453    c$$$     $              ,' >>> TRACK IMAGE >>> of'
454    c$$$     $              ,ibest
455    c$$$               goto 122         !image track found
456    c$$$            endif
457    c$$$         enddo
458    *     ---------------------------------------------------------------
459    *     take the candidate with the greatest number of matching couples
460    *     if more than one satisfies the requirement,
461    *     choose the one with more points and lower chi2
462    *     ---------------------------------------------------------------
463    *     count the number of matching couples
464             ncpmax = 0
465             iimage   = 0           !id of candidate with better matching
466           do i=1,ntracks           do i=1,ntracks
467              iimage=i              ncp=0
468              do ip=1,nplanes              do ip=1,nplanes
469                 if(     CP_STORE(nplanes-ip+1,icand).ne.                 if(CP_STORE(nplanes-ip+1,icand).ne.0)then
470       $              -1*CP_STORE(nplanes-ip+1,i) )iimage=0                    if(
471         $                 CP_STORE(nplanes-ip+1,i).ne.0
472         $                 .and.
473         $                 CP_STORE(nplanes-ip+1,icand).eq.
474         $                 -1*CP_STORE(nplanes-ip+1,i)
475         $                 )then
476                         ncp=ncp+1  !ok
477                      elseif(
478         $                    CP_STORE(nplanes-ip+1,i).ne.0
479         $                    .and.
480         $                    CP_STORE(nplanes-ip+1,icand).ne.
481         $                    -1*CP_STORE(nplanes-ip+1,i)
482         $                    )then
483                         ncp = 0
484                         goto 117   !it is not an image candidate
485                      else
486                        
487                      endif
488                   endif
489    c$$$               print*,' track ',i,' CP ',CP_STORE(nplanes-ip+1,i)
490    c$$$     $              ,CP_STORE(nplanes-ip+1,icand),ncp
491              enddo              enddo
492              if(  iimage.ne.0.and.   117        continue
493  c     $           RCHI2_STORE(i).le.CHI2MAX.and.              trackimage(i)=ncp   !number of matching couples
494  c     $           RCHI2_STORE(i).gt.0.and.              if(ncp>ncpmax)then
495       $           .true.)then                 ncpmax=ncp
496                 if(DEBUG)print*,'Track candidate ',iimage                 iimage=i
      $              ,' >>> TRACK IMAGE >>> of'  
      $              ,ibest  
                goto 122         !image track found  
497              endif              endif
498           enddo           enddo
499    *     check if there are more than one image candidates
500             ngood=0
501             do i=1,ntracks
502                if( ncpmax.ne.0.and.trackimage(i).eq.ncpmax )ngood=ngood+1
503             enddo
504    *     if there are, choose the best one
505             if(ngood.gt.1)then
506    *     -------------------------------------------------------
507    *     order track-candidates according to:
508    *     1st) decreasing n.points
509    *     2nd) increasing chi**2
510    *     -------------------------------------------------------
511                rchi2best=1000000000.
512                ndofbest=0            
513                do i=1,ntracks
514                   if( trackimage(i).eq.ncpmax )then
515                      ndof=0              
516                      do ii=1,nplanes    
517                         ndof=ndof        
518         $                    +int(xgood_store(ii,i))
519         $                    +int(ygood_store(ii,i))
520                      enddo              
521                      if(ndof.gt.ndofbest)then
522                         iimage=i
523                         rchi2best=RCHI2_STORE(i)
524                         ndofbest=ndof    
525                      elseif(ndof.eq.ndofbest)then
526                         if(RCHI2_STORE(i).lt.rchi2best.and.
527         $                    RCHI2_STORE(i).gt.0)then
528                            iimage=i
529                            rchi2best=RCHI2_STORE(i)
530                            ndofbest=ndof  
531                         endif            
532                      endif
533                   endif
534                enddo
535                
536             endif
537    
538             if(DEBUG.EQ.1)print*,'Track candidate ',iimage
539         $        ,' >>> TRACK IMAGE >>> of'
540         $        ,ibest
541    
542   122     continue   122     continue
543    
544    
545  *     --- and store the results --------------------------------  *     --- and store the results --------------------------------
546           ntrk = ntrk + 1                   !counter of found tracks           ntrk = ntrk + 1                   !counter of found tracks
547           if(.not.FIMAGE           if(.not.FIMAGE
# Line 312  c         print*,'++++++++++ iimage,fima Line 554  c         print*,'++++++++++ iimage,fima
554  c     $        ,iimage,fimage,ntrk,image(ntrk)  c     $        ,iimage,fimage,ntrk,image(ntrk)
555    
556           if(ntrk.eq.NTRKMAX)then           if(ntrk.eq.NTRKMAX)then
557              if(verbose)              if(verbose.eq.1)
558       $           print*,       $           print*,
559       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
560       $           'tracks exceeds vector dimension '       $           'tracks exceeds vector dimension '
# Line 348  cc            good2=.false. Line 590  cc            good2=.false.
590       $        rchi2best.le.CHI2MAX.and.       $        rchi2best.le.CHI2MAX.and.
591  c     $        rchi2best.lt.15..and.  c     $        rchi2best.lt.15..and.
592       $        .true.)then       $        .true.)then
593              if(DEBUG)then              if(DEBUG.EQ.1)then
594                 print*,'***** NEW SEARCH ****'                 print*,'***** NEW SEARCH ****'
595              endif              endif
596              goto 11111          !try new search              goto 11111          !try new search
# Line 362  c     $        rchi2best.lt.15..and. Line 604  c     $        rchi2best.lt.15..and.
604        end        end
605    
606    
   
   
 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$$$  
   
607                
608  ************************************************************  ************************************************************
609  ************************************************************  ************************************************************
# Line 558  c$$$ Line 628  c$$$
628  *     PFAy   - Position Finding Algorithm in y (COG2,ETA2,...)  *     PFAy   - Position Finding Algorithm in y (COG2,ETA2,...)
629  *     angx   - Projected angle in x  *     angx   - Projected angle in x
630  *     angy   - Projected angle in y  *     angy   - Projected angle in y
631    *     bfx    - x-component of magnetci field
632    *     bfy    - y-component of magnetic field
633  *  *
634  *     --------- COUPLES -------------------------------------------------------  *     --------- COUPLES -------------------------------------------------------
635  *     The couple defines a point in the space.  *     The couple defines a point in the space.
# Line 596  c$$$ Line 668  c$$$
668  *  *
669  *  *
670    
671        subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,angx,angy)        subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
672    
 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*****************************************************  
673                
674        include 'commontracker.f'        include 'commontracker.f'
       include 'calib.f'  
675        include 'level1.f'        include 'level1.f'
676          include 'calib.f'
677        include 'common_align.f'        include 'common_align.f'
678        include 'common_mech.f'        include 'common_mech.f'
679        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
       include 'common_resxy.f'  
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
680    
681        integer icx,icy           !X-Y cluster ID        integer icx,icy           !X-Y cluster ID
682        integer sensor        integer sensor
683        integer viewx,viewy        integer viewx,viewy
684        character*4 PFAx,PFAy     !PFA to be used        character*4 PFAx,PFAy     !PFA to be used
685        real angx,angy            !X-Y angle        real ax,ay                !X-Y geometric angle
686          real angx,angy            !X-Y effective angle
687          real bfx,bfy              !X-Y b-field components
688    
689        real stripx,stripy        real stripx,stripy
690    
691        double precision xrt,yrt,zrt        double precision xrt,yrt,zrt
692        double precision xrt_A,yrt_A,zrt_A        double precision xrt_A,yrt_A,zrt_A
693        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  
694                
695    
696        parameter (ndivx=30)        parameter (ndivx=30)
697    
698    
699    c$$$      print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy
700                
701        resxPAM = 0        resxPAM = 0
702        resyPAM = 0        resyPAM = 0
# Line 648  c      double precision xi_B,yi_B,zi_B Line 710  c      double precision xi_B,yi_B,zi_B
710        xPAM_B = 0.        xPAM_B = 0.
711        yPAM_B = 0.        yPAM_B = 0.
712        zPAM_B = 0.        zPAM_B = 0.
713    c      print*,'## xyz_PAM: ',icx,icy,sensor,PFAx,PFAy,angx,angy
714    
715          if(sensor.lt.1.or.sensor.gt.2)then
716             print*,'xyz_PAM   ***ERROR*** wrong input '
717             print*,'sensor ',sensor
718             icx=0
719             icy=0
720          endif
721    
722  *     -----------------  *     -----------------
723  *     CLUSTER X  *     CLUSTER X
724  *     -----------------  *     -----------------      
   
725        if(icx.ne.0)then        if(icx.ne.0)then
726           viewx = VIEW(icx)  
727           nldx = nld(MAXS(icx),VIEW(icx))           viewx   = VIEW(icx)
728           nplx = npl(VIEW(icx))           nldx    = nld(MAXS(icx),VIEW(icx))
729           resxPAM = RESXAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!           nplx    = npl(VIEW(icx))
730            c         resxPAM = RESXAV
731           stripx = float(MAXS(icx))           stripx  = float(MAXS(icx))
732           if(PFAx.eq.'COG1')then  !(1)  
733              stripx = stripx      !(1)           if(
734              resxPAM = resxPAM    !(1)       $        viewx.lt.1.or.        
735           elseif(PFAx.eq.'COG2')then       $        viewx.gt.12.or.
736              stripx = stripx + cog(2,icx)                   $        nldx.lt.1.or.
737              resxPAM = resxPAM*fbad_cog(2,icx)       $        nldx.gt.3.or.
738           elseif(PFAx.eq.'ETA2')then       $        stripx.lt.1.or.
739  c            cog2 = cog(2,icx)       $        stripx.gt.3072.or.
740  c            etacorr = pfaeta2(cog2,viewx,nldx,angx)                   $        .false.)then
741  c            stripx = stripx + etacorr              print*,'xyz_PAM   ***ERROR*** wrong input '
742              stripx = stripx + pfaeta2(icx,angx)            !(3)              print*,'icx ',icx,'view ',viewx,'nld ',nldx,'strip ',stripx
743              resxPAM = risx_eta2(angx)                       !   (4)              icx = 0
744              if(DEBUG.and.fbad_cog(2,icx).ne.1)              goto 10
      $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)  
             resxPAM = resxPAM*fbad_cog(2,icx)  
          elseif(PFAx.eq.'ETA3')then                         !(3)  
             stripx = stripx + pfaeta3(icx,angx)            !(3)  
             resxPAM = risx_eta3(angx)                       !   (4)  
             if(DEBUG.and.fbad_cog(3,icx).ne.1)              !(3)  
      $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)  
             resxPAM = resxPAM*fbad_cog(3,icx)               !(3)  
          elseif(PFAx.eq.'ETA4')then                         !(3)  
             stripx = stripx + pfaeta4(icx,angx)            !(3)  
             resxPAM = risx_eta4(angx)                       !   (4)  
             if(DEBUG.and.fbad_cog(4,icx).ne.1)              !(3)  
      $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)  
             resxPAM = resxPAM*fbad_cog(4,icx)               !(3)  
          elseif(PFAx.eq.'ETA')then                          !(3)  
             stripx = stripx + pfaeta(icx,angx)             !(3)  
             resxPAM = ris_eta(icx,angx)                     !   (4)  
             if(DEBUG.and.fbad_cog(2,icx).ne.1)              !(3)  
      $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)  
 c            resxPAM = resxPAM*fbad_cog(2,icx)               !(3)TEMPORANEO  
             resxPAM = resxPAM*fbad_eta(icx,angx)             !(3)(4)  
          elseif(PFAx.eq.'COG')then           !(2)  
             stripx = stripx + cog(0,icx)     !(2)          
             resxPAM = risx_cog(angx)                        !   (4)  
             resxPAM = resxPAM*fbad_cog(0,icx)!(2)  
          else  
             print*,'*** Non valid p.f.a. (x) --> ',PFAx  
745           endif           endif
746    
747        endif  *        --------------------------
748  c      if(icy.eq.0.and.icx.ne.0)  *        magnetic-field corrections
749  c     $     print*,PFAx,icx,angx,stripx,resxPAM,'***'  *        --------------------------
750                   stripx  = stripx +  fieldcorr(viewx,bfy)      
751             angx    = effectiveangle(ax,viewx,bfy)
752            
753             call applypfa(PFAx,icx,angx,corr,res)
754             stripx  = stripx + corr
755             resxPAM = res
756    
757     10   endif
758        
759  *     -----------------  *     -----------------
760  *     CLUSTER Y  *     CLUSTER Y
761  *     -----------------  *     -----------------
762    
763        if(icy.ne.0)then        if(icy.ne.0)then
764    
765           viewy = VIEW(icy)           viewy = VIEW(icy)
766           nldy = nld(MAXS(icy),VIEW(icy))           nldy = nld(MAXS(icy),VIEW(icy))
767           nply = npl(VIEW(icy))           nply = npl(VIEW(icy))
768           resyPAM = RESYAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!  c         resyPAM = RESYAV
769             stripy = float(MAXS(icy))
770    
771             if(
772         $        viewy.lt.1.or.        
773         $        viewy.gt.12.or.
774         $        nldy.lt.1.or.
775         $        nldy.gt.3.or.
776         $        stripy.lt.1.or.
777         $        stripy.gt.3072.or.
778         $        .false.)then
779                print*,'xyz_PAM   ***ERROR*** wrong input '
780                print*,'icy ',icy,'view ',viewy,'nld ',nldy,'strip ',stripy
781                icy = 0
782                goto 20
783             endif
784    
785           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
786              print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '              if(DEBUG.EQ.1) then
787       $           ,icx,icy                 print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '
788         $              ,icx,icy
789                endif
790              goto 100              goto 100
791           endif           endif
           
          stripy = float(MAXS(icy))  
          if(PFAy.eq.'COG1')then !(1)  
             stripy = stripy     !(1)  
             resyPAM = resyPAM   !(1)  
          elseif(PFAy.eq.'COG2')then  
             stripy = stripy + cog(2,icy)  
             resyPAM = resyPAM*fbad_cog(2,icy)  
          elseif(PFAy.eq.'ETA2')then  
 c            cog2 = cog(2,icy)  
 c            etacorr = pfaeta2(cog2,viewy,nldy,angy)  
 c            stripy = stripy + etacorr  
             stripy = stripy + pfaeta2(icy,angy)            !(3)  
             resyPAM = risy_eta2(angy)                       !   (4)  
             resyPAM = resyPAM*fbad_cog(2,icy)  
             if(DEBUG.and.fbad_cog(2,icy).ne.1)  
      $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)  
          elseif(PFAy.eq.'ETA3')then                         !(3)  
             stripy = stripy + pfaeta3(icy,angy)            !(3)  
             resyPAM = resyPAM*fbad_cog(3,icy)               !(3)  
             if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)  
      $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)  
          elseif(PFAy.eq.'ETA4')then                         !(3)  
             stripy = stripy + pfaeta4(icy,angy)            !(3)  
             resyPAM = resyPAM*fbad_cog(4,icy)               !(3)  
             if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)  
      $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)  
          elseif(PFAy.eq.'ETA')then                          !(3)  
             stripy = stripy + pfaeta(icy,angy)             !(3)  
             resyPAM = ris_eta(icy,angy)                     !   (4)  
 c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO  
             resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)  
             if(DEBUG.and.fbad_cog(2,icy).ne.1)              !(3)  
      $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)!(3)  
          elseif(PFAy.eq.'COG')then  
             stripy = stripy + cog(0,icy)              
             resyPAM = risy_cog(angy)                        !   (4)  
 c            resyPAM = ris_eta(icy,angy)                    !   (4)  
             resyPAM = resyPAM*fbad_cog(0,icy)  
          else  
             print*,'*** Non valid p.f.a. (x) --> ',PFAx  
          endif  
792    
793        endif  *        --------------------------
794    *        magnetic-field corrections
795    *        --------------------------
796             stripy  = stripy +  fieldcorr(viewy,bfx)      
797             angy    = effectiveangle(ay,viewy,bfx)
798            
799             call applypfa(PFAy,icy,angy,corr,res)
800             stripy  = stripy + corr
801             resyPAM = res
802    
803     20   endif
804    
805    c$$$      print*,'## stripx,stripy ',stripx,stripy
806    
         
807  c===========================================================  c===========================================================
808  C     COUPLE  C     COUPLE
809  C===========================================================  C===========================================================
# Line 778  c     (xi,yi,zi) = mechanical coordinate Line 814  c     (xi,yi,zi) = mechanical coordinate
814  c------------------------------------------------------------------------  c------------------------------------------------------------------------
815           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
816       $        .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...
817              print*,'xyz_PAM (couple):',              if(DEBUG.EQ.1) then
818       $          ' WARNING: false X strip: strip ',stripx                 print*,'xyz_PAM (couple):',
819         $              ' WARNING: false X strip: strip ',stripx
820                endif
821           endif           endif
822           xi = acoordsi(stripx,viewx)           xi = acoordsi(stripx,viewx)
823           yi = acoordsi(stripy,viewy)           yi = acoordsi(stripy,viewy)
824           zi = 0.           zi = 0.
825                    
   
826  c------------------------------------------------------------------------  c------------------------------------------------------------------------
827  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
828  c------------------------------------------------------------------------  c------------------------------------------------------------------------
# Line 871  c            print*,'X-singlet ',icx,npl Line 908  c            print*,'X-singlet ',icx,npl
908  c            if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...  c            if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...
909              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
910       $           .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...
911                 print*,'xyz_PAM (X-singlet):',                 if(DEBUG.EQ.1) then
912       $             ' WARNING: false X strip: strip ',stripx                    print*,'xyz_PAM (X-singlet):',
913         $                 ' WARNING: false X strip: strip ',stripx
914                   endif
915              endif              endif
916              xi   = acoordsi(stripx,viewx)              xi   = acoordsi(stripx,viewx)
917    
# Line 894  c            print*,'X-cl ',icx,stripx,' Line 933  c            print*,'X-cl ',icx,stripx,'
933  c            print*,yi_A,' <--> ',yi_B  c            print*,yi_A,' <--> ',yi_B
934    
935           else           else
936                if(DEBUG.EQ.1) then
937              print *,'routine xyz_PAM ---> not properly used !!!'                 print *,'routine xyz_PAM ---> not properly used !!!'
938              print *,'icx = ',icx                 print *,'icx = ',icx
939              print *,'icy = ',icy                 print *,'icy = ',icy
940                endif
941              goto 100              goto 100
942                            
943           endif           endif
# Line 962  c--------------------------------------- Line 1002  c---------------------------------------
1002  c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'  c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'
1003    
1004        else        else
1005                       if(DEBUG.EQ.1) then
1006           print *,'routine xyz_PAM ---> not properly used !!!'              print *,'routine xyz_PAM ---> not properly used !!!'
1007           print *,'icx = ',icx              print *,'icx = ',icx
1008           print *,'icy = ',icy              print *,'icy = ',icy
1009                         endif
1010        endif        endif
1011                    
1012    
1013    c      print*,'## xPAM,yPAM,zPAM       ',xPAM,yPAM,zPAM
1014    c      print*,'## xPAM_A,yPAM_A,zPAM_A ',xPAM_A,yPAM_A,zPAM_A
1015    c      print*,'## xPAM_B,yPAM_B,zPAM_B ',xPAM_B,yPAM_B,zPAM_B
1016    
1017   100  continue   100  continue
1018        end        end
1019    
1020    ************************************************************************
1021    *     Call xyz_PAM subroutine with default PFA and fill the mini2 common.
1022    *     (done to be called from c/c++)
1023    ************************************************************************
1024    
1025          subroutine xyzpam(ip,icx,icy,lad,sensor,ax,ay,bfx,bfy)
1026    
1027          include 'commontracker.f'
1028          include 'level1.f'
1029          include 'common_mini_2.f'
1030          include 'common_xyzPAM.f'
1031          include 'common_mech.f'
1032          include 'calib.f'
1033          
1034    *     flag to chose PFA
1035    c$$$      character*10 PFA
1036    c$$$      common/FINALPFA/PFA
1037    
1038          integer icx,icy           !X-Y cluster ID
1039          integer sensor
1040          character*4 PFAx,PFAy     !PFA to be used
1041          real ax,ay                !X-Y geometric angle
1042          real bfx,bfy              !X-Y b-field components
1043    
1044          ipx=0
1045          ipy=0      
1046          
1047    c$$$      PFAx = 'COG4'!PFA
1048    c$$$      PFAy = 'COG4'!PFA
1049    
1050    
1051          if(icx.gt.nclstr1.or.icy.gt.nclstr1)then
1052                print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1053         $           ,' does not exists (nclstr1=',nclstr1,')'
1054                icx = -1*icx
1055                icy = -1*icy
1056                return
1057            
1058          endif
1059          
1060          call idtoc(pfaid,PFAx)
1061          call idtoc(pfaid,PFAy)
1062    
1063    cc      print*,PFAx,PFAy
1064    
1065    c$$$      call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1066    
1067    c$$$      print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy
1068          
1069          if(icx.ne.0.and.icy.ne.0)then
1070    
1071             ipx=npl(VIEW(icx))
1072             ipy=npl(VIEW(icy))
1073    c$$$         if( (nplanes-ipx+1).ne.ip.or.(nplanes-ipy+1).ne.ip )
1074    c$$$     $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1075    c$$$     $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1076            
1077             if( (nplanes-ipx+1).ne.ip )then
1078                print*,'xyzpam: ***WARNING*** cluster ',icx
1079         $           ,' does not belong to plane: ',ip
1080                icx = -1*icx
1081                return
1082             endif
1083             if( (nplanes-ipy+1).ne.ip )then
1084                print*,'xyzpam: ***WARNING*** cluster ',icy
1085         $           ,' does not belong to plane: ',ip
1086                icy = -1*icy
1087                return
1088             endif
1089    
1090             call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1091    
1092             xgood(ip) = 1.
1093             ygood(ip) = 1.
1094             resx(ip)  = resxPAM
1095             resy(ip)  = resyPAM
1096    
1097             xm(ip) = xPAM
1098             ym(ip) = yPAM
1099             zm(ip) = zPAM
1100             xm_A(ip) = 0.
1101             ym_A(ip) = 0.
1102             xm_B(ip) = 0.
1103             ym_B(ip) = 0.
1104    
1105    c         zv(ip) = zPAM
1106    
1107          elseif(icx.eq.0.and.icy.ne.0)then
1108    
1109             ipy=npl(VIEW(icy))
1110    c$$$         if((nplanes-ipy+1).ne.ip)
1111    c$$$     $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1112    c$$$     $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1113             if( (nplanes-ipy+1).ne.ip )then
1114                print*,'xyzpam: ***WARNING*** cluster ',icy
1115         $           ,' does not belong to plane: ',ip
1116                icy = -1*icy
1117                return
1118             endif
1119    
1120             call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1121            
1122             xgood(ip) = 0.
1123             ygood(ip) = 1.
1124             resx(ip)  = 1000.
1125             resy(ip)  = resyPAM
1126    
1127             xm(ip) = -100.
1128             ym(ip) = -100.
1129             zm(ip) = (zPAM_A+zPAM_B)/2.
1130             xm_A(ip) = xPAM_A
1131             ym_A(ip) = yPAM_A
1132             xm_B(ip) = xPAM_B
1133             ym_B(ip) = yPAM_B
1134    
1135    c         zv(ip) = (zPAM_A+zPAM_B)/2.
1136            
1137          elseif(icx.ne.0.and.icy.eq.0)then
1138    
1139             ipx=npl(VIEW(icx))
1140    c$$$         if((nplanes-ipx+1).ne.ip)
1141    c$$$     $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1142    c$$$     $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1143    
1144             if( (nplanes-ipx+1).ne.ip )then
1145                print*,'xyzpam: ***WARNING*** cluster ',icx
1146         $           ,' does not belong to plane: ',ip
1147                icx = -1*icx
1148                return
1149             endif
1150    
1151             call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1152          
1153             xgood(ip) = 1.
1154             ygood(ip) = 0.
1155             resx(ip)  = resxPAM
1156             resy(ip)  = 1000.
1157    
1158             xm(ip) = -100.
1159             ym(ip) = -100.
1160             zm(ip) = (zPAM_A+zPAM_B)/2.
1161             xm_A(ip) = xPAM_A
1162             ym_A(ip) = yPAM_A
1163             xm_B(ip) = xPAM_B
1164             ym_B(ip) = yPAM_B
1165            
1166    c         zv(ip) = (zPAM_A+zPAM_B)/2.
1167    
1168          else
1169    
1170             il = 2
1171             if(lad.ne.0)il=lad
1172             is = 1
1173             if(sensor.ne.0)is=sensor
1174    c         print*,nplanes-ip+1,il,is
1175    
1176             xgood(ip) = 0.
1177             ygood(ip) = 0.
1178             resx(ip)  = 1000.
1179             resy(ip)  = 1000.
1180    
1181             xm(ip) = -100.
1182             ym(ip) = -100.          
1183             zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1184             xm_A(ip) = 0.
1185             ym_A(ip) = 0.
1186             xm_B(ip) = 0.
1187             ym_B(ip) = 0.
1188    
1189    c         zv(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1190    
1191          endif
1192    
1193          if(DEBUG.EQ.1)then
1194    c         print*,'----------------------------- track coord'
1195    22222    format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
1196             write(*,22222)ip,zm(ip),xm(ip),ym(ip)
1197         $        ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip)
1198         $        ,xgood(ip),ygood(ip),resx(ip),resy(ip)
1199    c$$$         print*,'-----------------------------'
1200          endif
1201          end
1202    
1203  ********************************************************************************  ********************************************************************************
1204  ********************************************************************************  ********************************************************************************
# Line 1047  c         print*,'A-(',xPAM_A,yPAM_A,') Line 1274  c         print*,'A-(',xPAM_A,yPAM_A,')
1274           endif                   endif        
1275    
1276           distance=           distance=
1277       $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2       $       ((xmi-XPP)**2+(ymi-YPP)**2)!QUIQUI
1278    cc     $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2
1279           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1280    
1281  c$$$         print*,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b  c$$$         print*,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b
# Line 1072  c$$$         print*,' resolution ',re Line 1300  c$$$         print*,' resolution ',re
1300  *     ----------------------  *     ----------------------
1301                    
1302           distance=           distance=
1303       $        ((xPAM-XPP)/resxPAM)**2       $       ((xPAM-XPP))**2 !QUIQUI
1304       $        +       $       +
1305       $        ((yPAM-YPP)/resyPAM)**2       $       ((yPAM-YPP))**2
1306    c$$$     $        ((xPAM-XPP)/resxPAM)**2
1307    c$$$     $        +
1308    c$$$     $        ((yPAM-YPP)/resyPAM)**2
1309           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1310    
1311  c$$$         print*,xPAM,yPAM,zPAM  c$$$         print*,xPAM,yPAM,zPAM
# Line 1083  c$$$         print*,' resolution ',resxP Line 1314  c$$$         print*,' resolution ',resxP
1314                    
1315        else        else
1316                    
1317           print*  c         print*
1318       $        ,' function distance_to ---> wrong usage!!!'  c     $        ,' function distance_to ---> wrong usage!!!'
1319           print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM  c         print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM
1320           print*,' xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b '  c         print*,' xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b '
1321       $        ,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b  c     $        ,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b
1322        endif          endif  
1323    
1324        distance_to = sngl(distance)        distance_to = sngl(distance)
# Line 1155  c--------------------------------------- Line 1386  c---------------------------------------
1386                 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)                 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
1387       $              .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...
1388  c     if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...  c     if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...
1389                    print*,'whichsensor: ',  c                  print*,'whichsensor: ',
1390       $                ' WARNING: false X strip: strip ',stripx  c     $                ' WARNING: false X strip: strip ',stripx
1391                 endif                 endif
1392                 xi = acoordsi(stripx,viewx)                 xi = acoordsi(stripx,viewx)
1393                 yi = acoordsi(stripy,viewy)                 yi = acoordsi(stripy,viewy)
# Line 1285  c     $              ,iv,xvv(iv),yvv(iv) Line 1516  c     $              ,iv,xvv(iv),yvv(iv)
1516  *     it returns the plane number  *     it returns the plane number
1517  *      *    
1518        include 'commontracker.f'        include 'commontracker.f'
1519          include 'level1.f'
1520  c      include 'common_analysis.f'  c      include 'common_analysis.f'
1521        include 'common_momanhough.f'        include 'common_momanhough.f'
1522                
# Line 1310  c      include 'common_analysis.f' Line 1542  c      include 'common_analysis.f'
1542        is_cp=0        is_cp=0
1543        if(id.lt.0)is_cp=1        if(id.lt.0)is_cp=1
1544        if(id.gt.0)is_cp=2        if(id.gt.0)is_cp=2
1545        if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'  c      if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'
1546    
1547        return        return
1548        end        end
# Line 1322  c      include 'common_analysis.f' Line 1554  c      include 'common_analysis.f'
1554  *     it returns the id number ON THE PLANE  *     it returns the id number ON THE PLANE
1555  *      *    
1556        include 'commontracker.f'        include 'commontracker.f'
1557          include 'level1.f'
1558  c      include 'common_analysis.f'  c      include 'common_analysis.f'
1559        include 'common_momanhough.f'        include 'common_momanhough.f'
1560                
# Line 1350  c      include 'common_analysis.f' Line 1583  c      include 'common_analysis.f'
1583  *     positive if sensor =2  *     positive if sensor =2
1584  *  *
1585        include 'commontracker.f'        include 'commontracker.f'
1586          include 'level1.f'
1587  c      include 'calib.f'  c      include 'calib.f'
1588  c      include 'level1.f'  c      include 'level1.f'
1589  c      include 'common_analysis.f'  c      include 'common_analysis.f'
# Line 1379  c      include 'common_analysis.f' Line 1613  c      include 'common_analysis.f'
1613  *************************************************************************  *************************************************************************
1614  *************************************************************************  *************************************************************************
1615  *************************************************************************  *************************************************************************
 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  
1616                
1617    
1618  ***************************************************  ***************************************************
# Line 1661  c$$$      end Line 1627  c$$$      end
1627        subroutine cl_to_couples(iflag)        subroutine cl_to_couples(iflag)
1628    
1629        include 'commontracker.f'        include 'commontracker.f'
1630          include 'level1.f'
1631        include 'common_momanhough.f'        include 'common_momanhough.f'
1632        include 'momanhough_init.f'  c      include 'momanhough_init.f'
1633        include 'calib.f'        include 'calib.f'
1634        include 'level1.f'  c      include 'level1.f'
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
1635    
1636  *     output flag  *     output flag
1637  *     --------------  *     --------------
# Line 1677  c      common/dbg/DEBUG Line 1641  c      common/dbg/DEBUG
1641        integer iflag        integer iflag
1642    
1643        integer badseed,badclx,badcly        integer badseed,badclx,badcly
1644        
1645          if(DEBUG.EQ.1)print*,'cl_to_couples:'
1646    
1647  *     init variables  *     init variables
1648        ncp_tot=0        ncp_tot=0
# Line 1692  c      common/dbg/DEBUG Line 1658  c      common/dbg/DEBUG
1658           ncls(ip)=0           ncls(ip)=0
1659        enddo        enddo
1660        do icl=1,nclstrmax_level2        do icl=1,nclstrmax_level2
1661           cl_single(icl)=1           cl_single(icl) = 1
1662           cl_good(icl)=0           cl_good(icl)   = 0
1663          enddo
1664          do iv=1,nviews
1665             ncl_view(iv)  = 0
1666             mask_view(iv) = 0      !all included
1667        enddo        enddo
1668                
1669  *     start association  *     count number of cluster per view
       ncouples=0  
       do icx=1,nclstr1          !loop on cluster (X)  
          if(mod(VIEW(icx),2).eq.1)goto 10  
           
 *     ----------------------------------------------------  
 *     cut on charge (X VIEW)  
 *     ----------------------------------------------------  
          if(dedx(icx).lt.dedx_x_min)then  
             cl_single(icx)=0  
             goto 10  
          endif  
 *     ----------------------------------------------------  
 *     cut on multiplicity (X VIEW)  
 *     ----------------------------------------------------  
          if(mult(icx).ge.mult_x_max)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  
          badclx=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  
             badclx=badclx*ibad  
          enddo  
 *     ----------------------------------------------------  
 *     >>> eliminato il taglio sulle BAD <<<  
 *     ----------------------------------------------------  
 c     if(badcl.eq.0)then  
 c     cl_single(icx)=0  
 c     goto 10  
 c     endif  
 *     ----------------------------------------------------  
           
          cl_good(icx)=1  
          nplx=npl(VIEW(icx))  
          nldx=nld(MAXS(icx),VIEW(icx))  
           
          do icy=1,nclstr1       !loop on cluster (Y)  
             if(mod(VIEW(icy),2).eq.0)goto 20  
               
 *     ----------------------------------------------------  
 *     cut on charge (Y VIEW)  
 *     ----------------------------------------------------  
             if(dedx(icy).lt.dedx_y_min)then  
                cl_single(icy)=0  
                goto 20  
             endif  
 *     ----------------------------------------------------  
 *     cut on multiplicity (X VIEW)  
 *     ----------------------------------------------------  
             if(mult(icy).ge.mult_y_max)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  
             badcly=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))  
                badcly=badcly*ibad  
             enddo  
 *     ----------------------------------------------------  
 *     >>> eliminato il taglio sulle BAD <<<  
 *     ----------------------------------------------------  
 c     if(badcl.eq.0)then  
 c     cl_single(icy)=0  
 c     goto 20  
 c     endif  
 *     ----------------------------------------------------  
               
             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  
 *     (modified to be applied only below saturation... obviously)  
   
 *     -------------------------------------------------------------  
 *     >>> eliminata (TEMPORANEAMENTE) la correlazione di carica <<<  
 *     -------------------------------------------------------------  
                if(  .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx)  
      $              .and.  
      $              .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx)  
      $              .and.  
      $              (badclx.eq.1.and.badcly.eq.1)  
      $              .and.  
      $              .true.)then  
   
                   ddd=(dedx(icy)  
      $                 -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))  
                   ddd=ddd/sqrt(kch(nplx,nldx)**2+1)  
   
 c                  cut = chcut * sch(nplx,nldx)  
   
                   sss=(kch(nplx,nldx)*dedx(icy)+dedx(icx)  
      $                 -kch(nplx,nldx)*cch(nplx,nldx))  
                   sss=sss/sqrt(kch(nplx,nldx)**2+1)  
                   cut = chcut * (16 + sss/50.)  
   
                   if(abs(ddd).gt.cut)then  
                      goto 20    !charge not consistent  
                   endif  
                endif  
                 
 *     ------------------> COUPLE <------------------  
 *     check to do not overflow vector dimentions  
 c$$$               if(ncp_plane(nplx).gt.ncouplemax)then  
 c$$$                  if(DEBUG)print*,  
 c$$$     $                    ' ** warning ** number of identified'//  
 c$$$     $                    ' couples on plane ',nplx,  
 c$$$     $                    ' exceeds vector dimention'//  
 c$$$     $                    ' ( ',ncouplemax,' )'  
 c$$$c     good2=.false.  
 c$$$c     goto 880   !fill ntp and go to next event  
 c$$$                  iflag=1  
 c$$$                  return  
 c$$$               endif  
                 
                if(ncp_plane(nplx).eq.ncouplemax)then  
                   if(verbose)print*,  
      $                 '** warning ** number of identified '//  
      $                 'couples on plane ',nplx,  
      $                 'exceeds vector dimention '  
      $                 ,'( ',ncouplemax,' )'  
 c     good2=.false.  
 c     goto 880   !fill ntp and go to next event                      
                   iflag=1  
                   return  
                endif  
                 
                ncp_plane(nplx) = ncp_plane(nplx) + 1  
                clx(nplx,ncp_plane(nplx))=icx  
                cly(nply,ncp_plane(nplx))=icy  
                cl_single(icx)=0  
                cl_single(icy)=0  
             endif                                
 *     ----------------------------------------------  
   
  20         continue  
          enddo                  !end loop on clusters(Y)  
           
  10      continue  
       enddo                     !end loop on clusters(X)  
         
         
1670        do icl=1,nclstr1        do icl=1,nclstr1
1671           if(cl_single(icl).eq.1)then           ncl_view(VIEW(icl)) = ncl_view(VIEW(icl)) + 1        
             ip=npl(VIEW(icl))  
             ncls(ip)=ncls(ip)+1  
             cls(ip,ncls(ip))=icl  
          endif  
1672        enddo        enddo
1673          *     mask views with too many clusters
1674                do iv=1,nviews
1675        if(DEBUG)then           if( ncl_view(iv).gt. nclusterlimit)then
1676           print*,'clusters  ',nclstr1  c            mask_view(iv) = 1
1677           print*,'good    ',(cl_good(i),i=1,nclstr1)              mask_view(iv) = mask_view(iv) + 2**0
1678           print*,'singles ',(cl_single(i),i=1,nclstr1)              if(DEBUG.EQ.1)
1679           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)       $        print*,' * WARNING * cl_to_couple: n.clusters > '
1680        endif       $        ,nclusterlimit,' on view ', iv,' --> masked!'
1681                   endif
       do ip=1,6  
          ncp_tot=ncp_tot+ncp_plane(ip)  
1682        enddo        enddo
 c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  
         
       if(ncp_tot.gt.ncp_max)then  
          if(verbose)print*,  
      $           '** warning ** number of identified '//  
      $           'couples exceeds upper limit for Hough tr. '  
      $           ,'( ',ncp_max,' )'              
 c            good2=.false.  
 c     goto 880       !fill ntp and go to next event  
          iflag=1  
          return  
       endif  
         
       return  
       end  
         
 ***************************************************  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 **************************************************  
       subroutine cl_to_couples_nocharge(iflag)  
1683    
       include 'commontracker.f'  
       include 'common_momanhough.f'  
       include 'momanhough_init.f'  
       include 'calib.f'  
       include 'level1.f'  
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
1684    
 *     output flag  
 *     --------------  
 *     0 = good event  
 *     1 = bad event  
 *     --------------  
       integer iflag  
   
       integer badseed,badcl  
   
 *     init variables  
       ncp_tot=0  
       do ip=1,nplanes  
          do ico=1,ncouplemax  
             clx(ip,ico)=0  
             cly(ip,ico)=0  
          enddo  
          ncp_plane(ip)=0  
          do icl=1,nclstrmax_level2  
             cls(ip,icl)=1  
          enddo  
          ncls(ip)=0  
       enddo  
       do icl=1,nclstrmax_level2  
          cl_single(icl)=1  
          cl_good(icl)=0  
       enddo  
         
1685  *     start association  *     start association
1686        ncouples=0        ncouples=0
1687        do icx=1,nclstr1          !loop on cluster (X)        do icx=1,nclstr1          !loop on cluster (X)
1688           if(mod(VIEW(icx),2).eq.1)goto 10           if(mod(VIEW(icx),2).eq.1)goto 10
1689                    
1690  *     ----------------------------------------------------  *     ----------------------------------------------------
1691    *     jump masked views (X VIEW)
1692    *     ----------------------------------------------------
1693             if( mask_view(VIEW(icx)).ne.0 ) goto 10
1694    *     ----------------------------------------------------
1695  *     cut on charge (X VIEW)  *     cut on charge (X VIEW)
1696           if(dedx(icx).lt.dedx_x_min)then  *     ----------------------------------------------------
1697             if(sgnl(icx).lt.dedx_x_min)then
1698              cl_single(icx)=0              cl_single(icx)=0
1699              goto 10              goto 10
1700           endif           endif
1701    *     ----------------------------------------------------
1702  *     cut BAD (X VIEW)              *     cut BAD (X VIEW)            
1703    *     ----------------------------------------------------
1704           badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))           badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))
1705           ifirst=INDSTART(icx)           ifirst=INDSTART(icx)
1706           if(icx.ne.nclstr1) then           if(icx.ne.nclstr1) then
# Line 1977  c      common/dbg/DEBUG Line 1708  c      common/dbg/DEBUG
1708           else           else
1709              ilast=TOTCLLENGTH              ilast=TOTCLLENGTH
1710           endif           endif
1711           badcl=badseed           badclx=badseed
1712           do igood=-ngoodstr,ngoodstr           do igood=-ngoodstr,ngoodstr
1713              ibad=1              ibad=1
1714              if((INDMAX(icx)+igood).gt.ifirst.and.              if((INDMAX(icx)+igood).gt.ifirst.and.
# Line 1987  c      common/dbg/DEBUG Line 1718  c      common/dbg/DEBUG
1718       $              nvk(MAXS(icx)+igood),       $              nvk(MAXS(icx)+igood),
1719       $              nst(MAXS(icx)+igood))       $              nst(MAXS(icx)+igood))
1720              endif              endif
1721              badcl=badcl*ibad              badclx=badclx*ibad
1722           enddo           enddo
1723           if(badcl.eq.0)then     !<<<<<<<<<<<<<< BAD cut  *     ----------------------------------------------------
1724              cl_single(icx)=0    !<<<<<<<<<<<<<< BAD cut  *     >>> eliminato il taglio sulle BAD <<<
1725              goto 10             !<<<<<<<<<<<<<< BAD cut  *     ----------------------------------------------------
1726           endif                  !<<<<<<<<<<<<<< BAD cut  c     if(badcl.eq.0)then
1727    c     cl_single(icx)=0
1728    c     goto 10
1729    c     endif
1730  *     ----------------------------------------------------  *     ----------------------------------------------------
1731                    
1732           cl_good(icx)=1           cl_good(icx)=1
# Line 2003  c      common/dbg/DEBUG Line 1737  c      common/dbg/DEBUG
1737              if(mod(VIEW(icy),2).eq.0)goto 20              if(mod(VIEW(icy),2).eq.0)goto 20
1738                            
1739  *     ----------------------------------------------------  *     ----------------------------------------------------
1740    *     jump masked views (Y VIEW)
1741    *     ----------------------------------------------------
1742                if( mask_view(VIEW(icy)).ne.0 ) goto 20
1743    
1744    *     ----------------------------------------------------
1745  *     cut on charge (Y VIEW)  *     cut on charge (Y VIEW)
1746              if(dedx(icy).lt.dedx_y_min)then  *     ----------------------------------------------------
1747                if(sgnl(icy).lt.dedx_y_min)then
1748                 cl_single(icy)=0                 cl_single(icy)=0
1749                 goto 20                 goto 20
1750              endif              endif
1751    *     ----------------------------------------------------
1752  *     cut BAD (Y VIEW)              *     cut BAD (Y VIEW)            
1753    *     ----------------------------------------------------
1754              badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))              badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))
1755              ifirst=INDSTART(icy)              ifirst=INDSTART(icy)
1756              if(icy.ne.nclstr1) then              if(icy.ne.nclstr1) then
# Line 2016  c      common/dbg/DEBUG Line 1758  c      common/dbg/DEBUG
1758              else              else
1759                 ilast=TOTCLLENGTH                 ilast=TOTCLLENGTH
1760              endif              endif
1761              badcl=badseed              badcly=badseed
1762              do igood=-ngoodstr,ngoodstr              do igood=-ngoodstr,ngoodstr
1763                 ibad=1                 ibad=1
1764                 if((INDMAX(icy)+igood).gt.ifirst.and.                 if((INDMAX(icy)+igood).gt.ifirst.and.
# Line 2025  c      common/dbg/DEBUG Line 1767  c      common/dbg/DEBUG
1767       $              ibad=BAD(VIEW(icy),       $              ibad=BAD(VIEW(icy),
1768       $              nvk(MAXS(icy)+igood),       $              nvk(MAXS(icy)+igood),
1769       $              nst(MAXS(icy)+igood))       $              nst(MAXS(icy)+igood))
1770                 badcl=badcl*ibad                 badcly=badcly*ibad
1771              enddo              enddo
             if(badcl.eq.0)then  !<<<<<<<<<<<<<< BAD cut  
                cl_single(icy)=0 !<<<<<<<<<<<<<< BAD cut  
                goto 20          !<<<<<<<<<<<<<< BAD cut  
             endif               !<<<<<<<<<<<<<< BAD cut  
1772  *     ----------------------------------------------------  *     ----------------------------------------------------
1773                *     >>> eliminato il taglio sulle BAD <<<
1774    *     ----------------------------------------------------
1775    c     if(badcl.eq.0)then
1776    c     cl_single(icy)=0
1777    c     goto 20
1778    c     endif
1779    *     ----------------------------------------------------
1780                            
1781              cl_good(icy)=1                                cl_good(icy)=1                  
1782              nply=npl(VIEW(icy))              nply=npl(VIEW(icy))
# Line 2043  c      common/dbg/DEBUG Line 1787  c      common/dbg/DEBUG
1787  *     ----------------------------------------------  *     ----------------------------------------------
1788  *     geometrical consistency (same plane and ladder)  *     geometrical consistency (same plane and ladder)
1789              if(nply.eq.nplx.and.nldy.eq.nldx)then              if(nply.eq.nplx.and.nldy.eq.nldx)then
1790  *     charge correlation  *     charge correlation
1791  *     ===========================================================  *     (modified to be applied only below saturation... obviously)
1792  *     this version of the subroutine is used for the calibration  
1793  *     thus charge-correlation selection is obviously removed                 if(  .not.(sgnl(icy).gt.chsaty.and.sgnl(icx).gt.chsatx)
1794  *     ===========================================================       $              .and.
1795  c$$$               ddd=(dedx(icy)       $              .not.(sgnl(icy).lt.chmipy.and.sgnl(icx).lt.chmipx)
1796  c$$$     $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))       $              .and.
1797  c$$$               ddd=ddd/sqrt(kch(nplx,nldx)**2+1)       $              (badclx.eq.1.and.badcly.eq.1)
1798  c$$$               cut=chcut*sch(nplx,nldx)       $              .and.
1799  c$$$               if(abs(ddd).gt.cut)goto 20 !charge not consistent       $              .true.)then
1800  *     ===========================================================  
1801                                    ddd=(sgnl(icy)
1802                       $                 -kch(nplx,nldx)*sgnl(icx)-cch(nplx,nldx))
1803  *     ------------------> COUPLE <------------------                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
1804  *     check to do not overflow vector dimentions  
1805  c$$$               if(ncp_plane(nplx).gt.ncouplemax)then  c                  cut = chcut * sch(nplx,nldx)
1806  c$$$                  if(DEBUG)print*,  
1807  c$$$     $                    ' ** warning ** number of identified'//                    sss=(kch(nplx,nldx)*sgnl(icy)+sgnl(icx)
1808  c$$$     $                    ' couples on plane ',nplx,       $                 -kch(nplx,nldx)*cch(nplx,nldx))
1809  c$$$     $                    ' exceeds vector dimention'//                    sss=sss/sqrt(kch(nplx,nldx)**2+1)
1810  c$$$     $                    ' ( ',ncouplemax,' )'                    cut = chcut * (16 + sss/50.)
1811  c$$$c     good2=.false.  
1812  c$$$c     goto 880   !fill ntp and go to next event                    if(abs(ddd).gt.cut)then
1813  c$$$                  iflag=1                       goto 20    !charge not consistent
1814  c$$$                  return                    endif
1815  c$$$               endif                 endif
1816                  
1817                 if(ncp_plane(nplx).eq.ncouplemax)then                 if(ncp_plane(nplx).gt.ncouplemax)then
1818                    if(verbose)print*,                    if(verbose.eq.1)print*,
1819       $                 '** warning ** number of identified '//       $                 '** warning ** number of identified '//
1820       $                 'couples on plane ',nplx,       $                 'couples on plane ',nplx,
1821       $                 'exceeds vector dimention '       $                 'exceeds vector dimention '
1822       $                 ,'( ',ncouplemax,' )'       $                 ,'( ',ncouplemax,' ) --> masked!'
1823  c     good2=.false.  c                  mask_view(nviewx(nplx)) = 2
1824  c     goto 880   !fill ntp and go to next event                      c                  mask_view(nviewy(nply)) = 2
1825                    iflag=1                    mask_view(nviewx(nplx))= mask_view(nviewx(nplx))+ 2**1
1826                    return                    mask_view(nviewy(nply))= mask_view(nviewy(nply))+ 2**1
1827                      goto 10
1828                 endif                 endif
1829                                
1830    *     ------------------> COUPLE <------------------
1831                 ncp_plane(nplx) = ncp_plane(nplx) + 1                 ncp_plane(nplx) = ncp_plane(nplx) + 1
1832                 clx(nplx,ncp_plane(nplx))=icx                 clx(nplx,ncp_plane(nplx))=icx
1833                 cly(nply,ncp_plane(nplx))=icy                 cly(nply,ncp_plane(nplx))=icy
1834                 cl_single(icx)=0                 cl_single(icx)=0
1835                 cl_single(icy)=0                 cl_single(icy)=0
             endif                                
1836  *     ----------------------------------------------  *     ----------------------------------------------
1837    
1838                endif                              
1839    
1840   20         continue   20         continue
1841           enddo                  !end loop on clusters(Y)           enddo                  !end loop on clusters(Y)
1842                    
# Line 2106  c     goto 880   !fill ntp and go to nex Line 1853  c     goto 880   !fill ntp and go to nex
1853        enddo        enddo
1854                
1855                
1856        if(DEBUG)then        if(DEBUG.EQ.1)then
1857           print*,'clusters  ',nclstr1           print*,'clusters  ',nclstr1
1858           print*,'good    ',(cl_good(i),i=1,nclstr1)           print*,'good    ',(cl_good(i),i=1,nclstr1)
1859           print*,'singles ',(cl_single(i),i=1,nclstr1)           print*,'singlets',(cl_single(i),i=1,nclstr1)
1860           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
1861        endif        endif
1862                
1863        do ip=1,6        do ip=1,6
1864           ncp_tot=ncp_tot+ncp_plane(ip)           ncp_tot = ncp_tot + ncp_plane(ip)
1865        enddo        enddo
 c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  
         
       if(ncp_tot.gt.ncp_max)then  
          if(verbose)print*,  
      $           '** warning ** number of identified '//  
      $           'couples exceeds upper limit for Hough tr. '  
      $           ,'( ',ncp_max,' )'              
 c            good2=.false.  
 c     goto 880       !fill ntp and go to next event  
          iflag=1  
          return  
       endif  
1866                
1867        return        return
1868        end        end
   
1869                
1870  ***************************************************  ***************************************************
1871  *                                                 *  *                                                 *
# Line 2143  c     goto 880       !fill ntp and go to Line 1877  c     goto 880       !fill ntp and go to
1877  **************************************************  **************************************************
1878    
1879        subroutine cp_to_doubtrip(iflag)        subroutine cp_to_doubtrip(iflag)
 c*****************************************************  
 c     02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
1880    
1881        include 'commontracker.f'        include 'commontracker.f'
1882          include 'level1.f'
1883        include 'common_momanhough.f'        include 'common_momanhough.f'
       include 'momanhough_init.f'  
1884        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
1885        include 'common_mini_2.f'        include 'common_mini_2.f'
1886        include 'calib.f'        include 'calib.f'
       include 'level1.f'  
1887    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
1888    
1889  *     output flag  *     output flag
1890  *     --------------  *     --------------
# Line 2188  c      double precision xm3,ym3,zm3 Line 1916  c      double precision xm3,ym3,zm3
1916        real xc,zc,radius        real xc,zc,radius
1917  *     -----------------------------  *     -----------------------------
1918    
1919          if(DEBUG.EQ.1)print*,'cp_to_doubtrip:'
1920    
1921    *     --------------------------------------------
1922    *     put a limit to the maximum number of couples
1923    *     per plane, in order to apply hough transform
1924    *     (couples recovered during track refinement)
1925    *     --------------------------------------------
1926          do ip=1,nplanes
1927             if(ncp_plane(ip).gt.ncouplelimit)then
1928    c            mask_view(nviewx(ip)) = 8
1929    c            mask_view(nviewy(ip)) = 8
1930                mask_view(nviewx(ip)) = mask_view(nviewx(ip)) + 2**7
1931                mask_view(nviewy(ip)) = mask_view(nviewy(ip)) + 2**7
1932             endif
1933          enddo
1934    
1935    
1936        ndblt=0                   !number of doublets        ndblt=0                   !number of doublets
1937        ntrpt=0                   !number of triplets        ntrpt=0                   !number of triplets
1938                
1939        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1
1940           do is1=1,2             !loop on sensors - COPPIA 1           if(  mask_view(nviewx(ip1)).ne.0 .or.
1941                     $        mask_view(nviewy(ip1)).ne.0 )goto 10 !skip plane
1942             do is1=1,2             !loop on sensors - COPPIA 1            
1943              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1
1944                 icx1=clx(ip1,icp1)                 icx1=clx(ip1,icp1)
1945                 icy1=cly(ip1,icp1)                 icy1=cly(ip1,icp1)
1946  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)
1947                 call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)  c               call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)
1948                   call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.)
1949                 xm1=xPAM                 xm1=xPAM
1950                 ym1=yPAM                 ym1=yPAM
1951                 zm1=zPAM                                   zm1=zPAM                  
1952  c     print*,'***',is1,xm1,ym1,zm1  c     print*,'***',is1,xm1,ym1,zm1
1953    
1954                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
1955                    do is2=1,2    !loop on sensors -ndblt COPPIA 2                    if(  mask_view(nviewx(ip2)).ne.0 .or.
1956                             $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane
1957                      do is2=1,2    !loop on sensors -ndblt COPPIA 2                    
1958                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
1959                          icx2=clx(ip2,icp2)                          icx2=clx(ip2,icp2)
1960                          icy2=cly(ip2,icp2)                          icy2=cly(ip2,icp2)
1961  c                        call xyz_PAM  c                        call xyz_PAM
1962  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)
1963    c                        call xyz_PAM
1964    c     $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)
1965                          call xyz_PAM                          call xyz_PAM
1966       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.)
1967                          xm2=xPAM                          xm2=xPAM
1968                          ym2=yPAM                          ym2=yPAM
1969                          zm2=zPAM                          zm2=zPAM
# Line 2224  c     $                       (icx2,icy2 Line 1973  c     $                       (icx2,icy2
1973  *     (2 couples needed)  *     (2 couples needed)
1974  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1975                          if(ndblt.eq.ndblt_max)then                          if(ndblt.eq.ndblt_max)then
1976                             if(verbose)print*,                             if(verbose.eq.1)print*,
1977       $                          '** warning ** number of identified '//       $                          '** warning ** number of identified '//
1978       $                          'doublets exceeds vector dimention '       $                          'doublets exceeds vector dimention '
1979       $                          ,'( ',ndblt_max,' )'       $                          ,'( ',ndblt_max,' )'
1980  c                           good2=.false.  c                           good2=.false.
1981  c                           goto 880 !fill ntp and go to next event  c                           goto 880 !fill ntp and go to next event
1982                               do iv=1,12
1983    c                              mask_view(iv) = 3
1984                                  mask_view(iv) = mask_view(iv)+ 2**2
1985                               enddo
1986                             iflag=1                             iflag=1
1987                             return                             return
1988                          endif                          endif
# Line 2263  c$$$ Line 2016  c$$$
2016  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2017    
2018    
2019                          if(ip2.eq.nplanes)goto 30 !no possible combination with 3 couples                          if(ip2.eq.nplanes)goto 31 !no possible combination with 3 couples
2020    
2021                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3
2022                               if(  mask_view(nviewx(ip3)).ne.0 .or.
2023         $                          mask_view(nviewy(ip3)).ne.0 )goto 30 !skip plane
2024                             do is3=1,2 !loop on sensors - COPPIA 3                             do is3=1,2 !loop on sensors - COPPIA 3
2025                                                                
2026                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3
# Line 2272  c$$$ Line 2028  c$$$
2028                                   icy3=cly(ip3,icp3)                                   icy3=cly(ip3,icp3)
2029  c                                 call xyz_PAM  c                                 call xyz_PAM
2030  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)
2031    c                                 call xyz_PAM
2032    c     $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)
2033                                   call xyz_PAM                                   call xyz_PAM
2034       $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)       $                                (icx3,icy3,is3,PFAdef,PFAdef
2035         $                                ,0.,0.,0.,0.)
2036                                   xm3=xPAM                                   xm3=xPAM
2037                                   ym3=yPAM                                   ym3=yPAM
2038                                   zm3=zPAM                                   zm3=zPAM
2039  *     find the circle passing through the three points  *     find the circle passing through the three points
2040    c$$$                                 call tricircle(3,xp,zp,angp,resp,chi
2041    c$$$     $                                ,xc,zc,radius,iflag)
2042                                     iflag = DEBUG
2043                                   call tricircle(3,xp,zp,angp,resp,chi                                   call tricircle(3,xp,zp,angp,resp,chi
2044       $                                ,xc,zc,radius,iflag)       $                                ,xc,zc,radius,iflag)
2045  c     print*,xc,zc,radius  c     print*,xc,zc,radius
# Line 2294  c     $                                 Line 2056  c     $                                
2056  *     (3 couples needed)  *     (3 couples needed)
2057  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2058                                   if(ntrpt.eq.ntrpt_max)then                                   if(ntrpt.eq.ntrpt_max)then
2059                                      if(verbose)print*,                                      if(verbose.eq.1)print*,
2060       $                     '** warning ** number of identified '//       $                     '** warning ** number of identified '//
2061       $                     'triplets exceeds vector dimention '       $                     'triplets exceeds vector dimention '
2062       $                    ,'( ',ntrpt_max,' )'       $                    ,'( ',ntrpt_max,' )'
2063  c                                    good2=.false.  c                                    good2=.false.
2064  c                                    goto 880 !fill ntp and go to next event  c                                    goto 880 !fill ntp and go to next event
2065                                        do iv=1,nviews
2066    c                                       mask_view(iv) = 4
2067                                           mask_view(iv)=mask_view(iv)+ 2**3
2068                                        enddo
2069                                      iflag=1                                      iflag=1
2070                                      return                                      return
2071                                   endif                                   endif
# Line 2338  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp Line 2104  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp
2104                                endif                                endif
2105                             enddo !end loop on COPPIA 3                             enddo !end loop on COPPIA 3
2106                          enddo   !end loop on sensors - COPPIA 3                          enddo   !end loop on sensors - COPPIA 3
2107     30                     continue
2108                       enddo      !end loop on planes  - COPPIA 3                       enddo      !end loop on planes  - COPPIA 3
2109   30                  continue   31                  continue
2110                                            
2111   1                enddo         !end loop on COPPIA 2   1                enddo         !end loop on COPPIA 2
2112                 enddo            !end loop on sensors - COPPIA 2                 enddo            !end loop on sensors - COPPIA 2
2113     20            continue
2114              enddo               !end loop on planes  - COPPIA 2              enddo               !end loop on planes  - COPPIA 2
2115                            
2116           enddo                  !end loop on COPPIA1           enddo                  !end loop on COPPIA1
2117        enddo                     !end loop on sensors - COPPIA 1        enddo                     !end loop on sensors - COPPIA 1
2118     10   continue
2119        enddo                     !end loop on planes  - COPPIA 1        enddo                     !end loop on planes  - COPPIA 1
2120                
2121        if(DEBUG)then        if(DEBUG.EQ.1)then
2122           print*,'--- doublets ',ndblt           print*,'--- doublets ',ndblt
2123           print*,'--- triplets ',ntrpt           print*,'--- triplets ',ntrpt
2124        endif        endif
# Line 2374  c     goto 880               !ntp fill Line 2143  c     goto 880               !ntp fill
2143        subroutine doub_to_YZcloud(iflag)        subroutine doub_to_YZcloud(iflag)
2144    
2145        include 'commontracker.f'        include 'commontracker.f'
2146          include 'level1.f'
2147        include 'common_momanhough.f'        include 'common_momanhough.f'
2148        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2149    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2150    
2151  *     output flag  *     output flag
2152  *     --------------  *     --------------
# Line 2397  c      common/dbg/DEBUG Line 2165  c      common/dbg/DEBUG
2165        integer cp_useds1(ncouplemaxtot) ! sensor 1        integer cp_useds1(ncouplemaxtot) ! sensor 1
2166        integer cp_useds2(ncouplemaxtot) ! sensor 2        integer cp_useds2(ncouplemaxtot) ! sensor 2
2167    
2168          if(DEBUG.EQ.1)print*,'doub_to_YZcloud:'
2169    
2170  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2171  *     classification of DOUBLETS  *     classification of DOUBLETS
# Line 2410  c      common/dbg/DEBUG Line 2179  c      common/dbg/DEBUG
2179        distance=0        distance=0
2180        nclouds_yz=0              !number of clouds        nclouds_yz=0              !number of clouds
2181        npt_tot=0        npt_tot=0
2182          nloop=0                  
2183     90   continue                  
2184        do idb1=1,ndblt           !loop (1) on DOUBLETS        do idb1=1,ndblt           !loop (1) on DOUBLETS
2185           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
2186                            
# Line 2513  c     print*,'*   idbref,idb2 ',idbref,i Line 2284  c     print*,'*   idbref,idb2 ',idbref,i
2284              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
2285           enddo           enddo
2286  c     print*,'>>>> ',ncpused,npt,nplused  c     print*,'>>>> ',ncpused,npt,nplused
2287           if(ncpused.lt.ncpyz_min)goto 2228 !next doublet  c         if(ncpused.lt.ncpyz_min)goto 2228 !next doublet
2288           if(npt.lt.nptyz_min)goto 2228 !next doublet           if(npt.lt.nptyz_min)goto 2228 !next doublet
2289           if(nplused.lt.nplyz_min)goto 2228 !next doublet           if(nplused.lt.nplyz_min)goto 2228 !next doublet
2290                    
# Line 2521  c     print*,'>>>> ',ncpused,npt,nplused Line 2292  c     print*,'>>>> ',ncpused,npt,nplused
2292  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2293    
2294           if(nclouds_yz.ge.ncloyz_max)then           if(nclouds_yz.ge.ncloyz_max)then
2295              if(verbose)print*,              if(verbose.eq.1)print*,
2296       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2297       $           'YZ clouds exceeds vector dimention '       $           'YZ clouds exceeds vector dimention '
2298       $           ,'( ',ncloyz_max,' )'       $           ,'( ',ncloyz_max,' )'
2299  c               good2=.false.  c               good2=.false.
2300  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2301                do iv=1,nviews
2302    c               mask_view(iv) = 5
2303                   mask_view(iv) = mask_view(iv) + 2**4
2304                enddo
2305              iflag=1              iflag=1
2306              return              return
2307           endif           endif
# Line 2545  c     ptcloud_yz_nt(nclouds_yz)=npt Line 2320  c     ptcloud_yz_nt(nclouds_yz)=npt
2320  c     print*,'>> ',ipt,db_all(ipt)  c     print*,'>> ',ipt,db_all(ipt)
2321           enddo             enddo  
2322           npt_tot=npt_tot+npt           npt_tot=npt_tot+npt
2323           if(DEBUG)then           if(DEBUG.EQ.1)then
2324              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
2325              print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'              print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'
2326              print*,'- alfayz1 ',alfayz1_av(nclouds_yz)              print*,'- alfayz1  ',alfayz1_av(nclouds_yz)
2327              print*,'- alfayz2 ',alfayz2_av(nclouds_yz)              print*,'- alfayz2  ',alfayz2_av(nclouds_yz)
2328              print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)              print*,'cp_useds1  ',(cp_useds1(icp),icp=1,ncp_tot)
2329              print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)              print*,'cp_useds2  ',(cp_useds2(icp),icp=1,ncp_tot)
2330              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)              print*,'cpcloud_yz '
2331         $           ,(cpcloud_yz(nclouds_yz,icp),icp=1,ncp_tot)
2332                print*,'hit_plane  ',(hit_plane(ip),ip=1,nplanes)
2333  c$$$            print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = '  c$$$            print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = '
2334  c$$$     $           ,ptcloud_yz(nclouds_yz)  c$$$     $           ,ptcloud_yz(nclouds_yz)
2335  c$$$            print*,'nt-uple: db_cloud(...) = '  c$$$            print*,'nt-uple: db_cloud(...) = '
# Line 2564  c$$$     $           ,(db_cloud(iii),iii Line 2341  c$$$     $           ,(db_cloud(iii),iii
2341        enddo                     !end loop (1) on DOUBLETS        enddo                     !end loop (1) on DOUBLETS
2342                
2343                
2344        if(DEBUG)then        if(nloop.lt.nstepy)then      
2345            cutdistyz = cutdistyz+cutystep
2346            nloop     = nloop+1          
2347            goto 90                
2348          endif                    
2349          
2350          if(DEBUG.EQ.1)then
2351           print*,'---------------------- '           print*,'---------------------- '
2352           print*,'Y-Z total clouds ',nclouds_yz           print*,'Y-Z total clouds ',nclouds_yz
2353           print*,' '           print*,' '
# Line 2590  c$$$     $           ,(db_cloud(iii),iii Line 2373  c$$$     $           ,(db_cloud(iii),iii
2373        subroutine trip_to_XZcloud(iflag)        subroutine trip_to_XZcloud(iflag)
2374    
2375        include 'commontracker.f'        include 'commontracker.f'
2376          include 'level1.f'
2377        include 'common_momanhough.f'        include 'common_momanhough.f'
2378        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2379    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2380    
2381  *     output flag  *     output flag
2382  *     --------------  *     --------------
# Line 2614  c      common/dbg/DEBUG Line 2396  c      common/dbg/DEBUG
2396        integer cp_useds1(ncouplemaxtot) ! sensor 1        integer cp_useds1(ncouplemaxtot) ! sensor 1
2397        integer cp_useds2(ncouplemaxtot) ! sensor 2        integer cp_useds2(ncouplemaxtot) ! sensor 2
2398    
2399          if(DEBUG.EQ.1)print*,'trip_to_XZcloud:'
2400    
2401  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2402  *     classification of TRIPLETS  *     classification of TRIPLETS
2403  *     according to distance in parameter space  *     according to distance in parameter space
# Line 2625  c      common/dbg/DEBUG Line 2409  c      common/dbg/DEBUG
2409        distance=0        distance=0
2410        nclouds_xz=0              !number of clouds                nclouds_xz=0              !number of clouds        
2411        npt_tot=0                 !total number of selected triplets        npt_tot=0                 !total number of selected triplets
2412          nloop=0                  
2413     91   continue                  
2414        do itr1=1,ntrpt           !loop (1) on TRIPLETS        do itr1=1,ntrpt           !loop (1) on TRIPLETS
2415           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
2416  c     print*,'--------------'  c     print*,'--------------'
# Line 2669  c         tr_temp(1)=itr1 Line 2455  c         tr_temp(1)=itr1
2455       $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2                     $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2              
2456                 distance = sqrt(distance)                 distance = sqrt(distance)
2457                                
2458                 if(distance.lt.cutdistxz)then  *     ------------------------------------------------------------------------
2459    *     FORCE INCLUSION OF TRIPLETS COMPOSED BY SAME COUPLES, IGNORING THE IMAGE
2460    *     ------------------------------------------------------------------------
2461    *     (added in august 2007)
2462                   istrimage=0
2463                   if(
2464         $              abs(cpxz1(itrref)).eq.abs(cpxz1(itr2)).and.
2465         $              abs(cpxz2(itrref)).eq.abs(cpxz2(itr2)).and.
2466         $              abs(cpxz3(itrref)).eq.abs(cpxz3(itr2)).and.
2467         $              .true.)istrimage=1
2468    
2469                   if(distance.lt.cutdistxz.or.istrimage.eq.1)then
2470  c     print*,idb1,idb2,distance,' cloud ',nclouds_yz  c     print*,idb1,idb2,distance,' cloud ',nclouds_yz
2471                    if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1                    if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1
2472                    if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1                    if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1
# Line 2726  c     print*,'check cp_used' Line 2523  c     print*,'check cp_used'
2523           do ip=1,nplanes           do ip=1,nplanes
2524              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
2525           enddo           enddo
2526           if(ncpused.lt.ncpxz_min)goto 22288 !next triplet  c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet
2527           if(npt.lt.nptxz_min)goto 22288     !next triplet           if(npt.lt.nptxz_min)goto 22288     !next triplet
2528           if(nplused.lt.nplxz_min)goto 22288 !next doublet           if(nplused.lt.nplxz_min)goto 22288 !next triplet
2529                    
2530  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2531  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2532           if(nclouds_xz.ge.ncloxz_max)then           if(nclouds_xz.ge.ncloxz_max)then
2533              if(verbose)print*,              if(verbose.eq.1)print*,
2534       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2535       $           'XZ clouds exceeds vector dimention '       $           'XZ clouds exceeds vector dimention '
2536       $           ,'( ',ncloxz_max,' )'       $           ,'( ',ncloxz_max,' )'
2537  c     good2=.false.  c     good2=.false.
2538  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2539                do iv=1,nviews
2540    c               mask_view(iv) = 6
2541                   mask_view(iv) =  mask_view(iv) + 2**5
2542                enddo
2543              iflag=1              iflag=1
2544              return              return
2545           endif           endif
# Line 2756  c     goto 880         !fill ntp and go Line 2557  c     goto 880         !fill ntp and go
2557           enddo           enddo
2558           npt_tot=npt_tot+npt           npt_tot=npt_tot+npt
2559                    
2560           if(DEBUG)then           if(DEBUG.EQ.1)then
2561              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
2562              print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'                            print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'              
2563              print*,'- alfaxz1 ',alfaxz1_av(nclouds_xz)              print*,'- alfaxz1  ',alfaxz1_av(nclouds_xz)
2564              print*,'- alfaxz2 ',alfaxz2_av(nclouds_xz)              print*,'- alfaxz2  ',alfaxz2_av(nclouds_xz)
2565              print*,'- alfaxz3 ',alfaxz3_av(nclouds_xz)              print*,'- alfaxz3  ',alfaxz3_av(nclouds_xz)
2566              print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)              print*,'cp_useds1  ',(cp_useds1(icp),icp=1,ncp_tot)
2567              print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)              print*,'cp_useds2  ',(cp_useds2(icp),icp=1,ncp_tot)
2568                print*,'cpcloud_xz '
2569         $           ,(cpcloud_xz(nclouds_xz,icp),icp=1,ncp_tot)
2570              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)
2571  c$$$            print*,'nt-uple: ptcloud_xz(',nclouds_xz,') = '  c$$$            print*,'nt-uple: ptcloud_xz(',nclouds_xz,') = '
2572  c$$$     $           ,ptcloud_xz(nclouds_xz)  c$$$     $           ,ptcloud_xz(nclouds_xz)
# Line 2774  c$$$     $           ,(tr_cloud(iii),iii Line 2577  c$$$     $           ,(tr_cloud(iii),iii
2577  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2578  22288    continue  22288    continue
2579        enddo                     !end loop (1) on DOUBLETS        enddo                     !end loop (1) on DOUBLETS
2580          
2581        if(DEBUG)then         if(nloop.lt.nstepx)then      
2582             cutdistxz=cutdistxz+cutxstep
2583             nloop=nloop+1          
2584             goto 91                
2585           endif                    
2586          
2587          if(DEBUG.EQ.1)then
2588           print*,'---------------------- '           print*,'---------------------- '
2589           print*,'X-Z total clouds ',nclouds_xz           print*,'X-Z total clouds ',nclouds_xz
2590           print*,' '           print*,' '
# Line 2796  c$$$     $           ,(tr_cloud(iii),iii Line 2605  c$$$     $           ,(tr_cloud(iii),iii
2605  **************************************************  **************************************************
2606    
2607        subroutine clouds_to_ctrack(iflag)        subroutine clouds_to_ctrack(iflag)
 c*****************************************************  
 c     02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
2608    
2609        include 'commontracker.f'        include 'commontracker.f'
2610          include 'level1.f'
2611        include 'common_momanhough.f'        include 'common_momanhough.f'
2612        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2613        include 'common_mini_2.f'        include 'common_mini_2.f'
2614        include 'common_mech.f'        include 'common_mech.f'
       include 'momanhough_init.f'  
2615    
2616  c      logical DEBUG  
 c      common/dbg/DEBUG  
2617    
2618  *     output flag  *     output flag
2619  *     --------------  *     --------------
# Line 2824  c      common/dbg/DEBUG Line 2629  c      common/dbg/DEBUG
2629  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2630  *     list of matching couples in the combination  *     list of matching couples in the combination
2631  *     between a XZ and YZ cloud  *     between a XZ and YZ cloud
2632        integer cp_match(nplanes,ncouplemax)        integer cp_match(nplanes,2*ncouplemax)
2633        integer ncp_match(nplanes)        integer ncp_match(nplanes)
2634  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2635        integer hit_plane(nplanes)        integer hit_plane(nplanes)
2636  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2637  *     variables for track fitting  *     variables for track fitting
2638        double precision AL_INI(5)        double precision AL_INI(5)
       double precision tath  
2639  *     -----------------------------------------------------------  *     -----------------------------------------------------------
 c      real fitz(nplanes)        !z coordinates of the planes in cm  
2640    
2641          if(DEBUG.EQ.1)print*,'clouds_to_ctrack:'
2642    
2643    
2644        ntracks=0                 !counter of track candidates        ntracks=0                 !counter of track candidates
# Line 2856  c      real fitz(nplanes)        !z coor Line 2660  c      real fitz(nplanes)        !z coor
2660                 enddo                 enddo
2661              enddo              enddo
2662              ncp_ok=0              ncp_ok=0
2663              do icp=1,ncp_tot    !loop on couples              do icp=1,ncp_tot    !loop over couples
2664  *     get info on  *     get info on
2665                 cpintersec(icp)=min(                 cpintersec(icp)=min(
2666       $              cpcloud_yz(iyz,icp),       $              cpcloud_yz(iyz,icp),
# Line 2865  c      real fitz(nplanes)        !z coor Line 2669  c      real fitz(nplanes)        !z coor
2669       $    (cpcloud_yz(iyz,icp).eq.1.and.cpcloud_xz(ixz,icp).eq.2).or.       $    (cpcloud_yz(iyz,icp).eq.1.and.cpcloud_xz(ixz,icp).eq.2).or.
2670       $    (cpcloud_yz(iyz,icp).eq.2.and.cpcloud_xz(ixz,icp).eq.1).or.       $    (cpcloud_yz(iyz,icp).eq.2.and.cpcloud_xz(ixz,icp).eq.1).or.
2671       $              .false.)cpintersec(icp)=0       $              .false.)cpintersec(icp)=0
2672    *     cpintersec is >0 if yz and xz clouds contain the same image of couple icp
2673                 if(cpintersec(icp).ne.0)then                 if(cpintersec(icp).ne.0)then
2674                    ncp_ok=ncp_ok+1                      ncp_ok=ncp_ok+1  
2675                                        
# Line 2897  c      real fitz(nplanes)        !z coor Line 2702  c      real fitz(nplanes)        !z coor
2702                 nplused=nplused+ hit_plane(ip)                 nplused=nplused+ hit_plane(ip)
2703              enddo              enddo
2704                            
             if(nplused.lt.nplxz_min)goto 888 !next doublet  
             if(ncp_ok.lt.ncpok)goto 888 !next cloud  
2705                            
2706              if(DEBUG)then              if(DEBUG.EQ.1)then
2707                 print*,'Combination ',iyz,ixz                 print*,'Combination ',iyz,ixz
2708       $              ,' db ',ptcloud_yz(iyz)       $              ,' db ',ptcloud_yz(iyz)
2709       $              ,' tr ',ptcloud_xz(ixz)       $              ,' tr ',ptcloud_xz(ixz)
2710       $              ,'  -----> # matching couples ',ncp_ok       $              ,'  -----> # matching couples ',ncp_ok
2711              endif              endif
2712    
2713    c            if(nplused.lt.nplxz_min)goto 888 !next combination
2714                if(nplused.lt.nplyz_min)goto 888 !next combination
2715                if(ncp_ok.lt.ncpok)goto 888 !next combination
2716    
2717  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
2718  c$$$  print*,'Configurazione cluster XZ'  c$$$  print*,'Configurazione cluster XZ'
2719  c$$$  print*,'1 -- ',(clx(1,i),i=1,ncp_plane(1))  c$$$  print*,'1 -- ',(clx(1,i),i=1,ncp_plane(1))
# Line 2924  c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_p Line 2732  c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_p
2732  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
2733                            
2734  *     -------> INITIAL GUESS <-------  *     -------> INITIAL GUESS <-------
2735              AL_INI(1)=dreal(alfaxz1_av(ixz))  cccc       SBAGLIATO
2736              AL_INI(2)=dreal(alfayz1_av(iyz))  c$$$            AL_INI(1) = dreal(alfaxz1_av(ixz))
2737              AL_INI(4)=datan(dreal(alfayz2_av(iyz))  c$$$            AL_INI(2) = dreal(alfayz1_av(iyz))
2738       $           /dreal(alfaxz2_av(ixz)))  c$$$            AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
2739              tath=-dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))  c$$$     $           /dreal(alfaxz2_av(ixz)))
2740              AL_INI(3)=tath/sqrt(1+tath**2)  c$$$            tath      = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
2741              AL_INI(5)=(1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.  c$$$            AL_INI(3) = tath/sqrt(1+tath**2)
2742                c$$$            AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
2743  c     print*,'*******',AL_INI(5)  cccc       GIUSTO (ma si sua guess())
2744              if(AL_INI(5).gt.defmax)goto 888 !next cloud  c$$$            AL_INI(1) = dreal(alfaxz1_av(ixz))
2745                c$$$            AL_INI(2) = dreal(alfayz1_av(iyz))
2746  c     print*,'alfaxz2, alfayz2 '  c$$$            tath      = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
2747  c     $              ,alfaxz2_av(ixz),alfayz2_av(iyz)  c$$$            AL_INI(3) = tath/sqrt(1+tath**2)
2748                c$$$            IF(alfaxz2_av(ixz).NE.0)THEN
2749  *     -------> INITIAL GUESS <-------  c$$$            AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
2750  c     print*,'AL_INI ',(al_ini(i),i=1,5)  c$$$     $           /dreal(alfaxz2_av(ixz)))
2751                c$$$            ELSE
2752              if(DEBUG)then  c$$$               AL_INI(4) = acos(-1.)/2
2753    c$$$               IF(alfayz2_av(iyz).LT.0)AL_INI(4) = AL_INI(4)+acos(-1.)
2754    c$$$            ENDIF
2755    c$$$            IF(alfaxz2_av(ixz).LT.0)AL_INI(4)= acos(-1.)+ AL_INI(4)
2756    c$$$            AL_INI(4) = -acos(-1.) + AL_INI(4) !from incidence direction to tracking rs
2757    c$$$            
2758    c$$$            AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
2759    c$$$            
2760    c$$$            if(AL_INI(5).gt.defmax)goto 888 !next cloud
2761                            
2762                if(DEBUG.EQ.1)then
2763                   print*,'track candidate', ntracks+1
2764                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))
2765                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))
2766                 print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4))                 print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4))
# Line 2974  c     print*,'AL_INI ',(al_ini(i),i=1,5) Line 2793  c     print*,'AL_INI ',(al_ini(i),i=1,5)
2793                                hit_plane(6)=icp6                                hit_plane(6)=icp6
2794                                if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6                                if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6
2795                                                                
2796                                  *                             ---------------------------------------
2797    *                             check if this group of couples has been
2798    *                             already fitted    
2799    *                             ---------------------------------------
2800                                  do ica=1,ntracks
2801                                     isthesame=1
2802                                     do ip=1,NPLANES
2803                                        if(hit_plane(ip).ne.0)then
2804                                           if(  CP_STORE(nplanes-ip+1,ica)
2805         $                                      .ne.
2806         $                                      cp_match(ip,hit_plane(ip)) )
2807         $                                      isthesame=0
2808                                        else
2809                                           if(  CP_STORE(nplanes-ip+1,ica)
2810         $                                      .ne.
2811         $                                      0 )
2812         $                                      isthesame=0
2813                                        endif
2814                                     enddo
2815                                     if(isthesame.eq.1)then
2816                                        if(DEBUG.eq.1)
2817         $                                   print*,'(already fitted)'
2818                                        goto 666 !jump to next combination
2819                                     endif
2820                                  enddo
2821    
2822                                call track_init !init TRACK common                                call track_init !init TRACK common
2823    
2824                                do ip=1,nplanes !loop on planes                                do ip=1,nplanes !loop on planes (bottom to top)
2825                                   if(hit_plane(ip).ne.0)then                                   if(hit_plane(ip).ne.0)then
2826                                      id=cp_match(ip,hit_plane(ip))                                      id=cp_match(ip,hit_plane(ip))
2827                                      is=is_cp(id)                                      is=is_cp(id)
# Line 2991  c     print*,'AL_INI ',(al_ini(i),i=1,5) Line 2835  c     print*,'AL_INI ',(al_ini(i),i=1,5)
2835  *                                   *************************  *                                   *************************
2836  c                                    call xyz_PAM(icx,icy,is,  c                                    call xyz_PAM(icx,icy,is,
2837  c     $                                   'COG2','COG2',0.,0.)  c     $                                   'COG2','COG2',0.,0.)
2838    c                                    call xyz_PAM(icx,icy,is, !(1)
2839    c     $                                   PFAdef,PFAdef,0.,0.) !(1)
2840                                      call xyz_PAM(icx,icy,is, !(1)                                      call xyz_PAM(icx,icy,is, !(1)
2841       $                                   PFAdef,PFAdef,0.,0.) !(1)       $                                   PFAdef,PFAdef,0.,0.,0.,0.)
2842  *                                   *************************  *                                   *************************
2843  *                                   -----------------------------  *                                   -----------------------------
2844                                      xgood(nplanes-ip+1)=1.                                      xgood(nplanes-ip+1)=1.
# Line 3008  c     $                                 Line 2854  c     $                                
2854  *     **********************************************************  *     **********************************************************
2855  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
2856  *     **********************************************************  *     **********************************************************
2857    cccc  scommentare se si usa al_ini della nuvola
2858    c$$$                              do i=1,5
2859    c$$$                                 AL(i)=AL_INI(i)
2860    c$$$                              enddo
2861                                  call guess()
2862                                do i=1,5                                do i=1,5
2863                                   AL(i)=AL_INI(i)                                   AL_INI(i)=AL(i)
2864                                enddo                                enddo
2865                                ifail=0 !error flag in chi^2 computation                                ifail=0 !error flag in chi^2 computation
2866                                jstep=0 !number of  minimization steps                                jstep=0 !number of  minimization steps
2867                                iprint=0                                iprint=0
2868                                if(DEBUG)iprint=1  c                              if(DEBUG.EQ.1)iprint=1
2869                                  if(DEBUG.EQ.1)iprint=2
2870                                call mini2(jstep,ifail,iprint)                                call mini2(jstep,ifail,iprint)
2871                                if(ifail.ne.0) then                                if(ifail.ne.0) then
2872                                   if(DEBUG)then                                   if(DEBUG.EQ.1)then
2873                                      print *,                                      print *,
2874       $                              '*** MINIMIZATION FAILURE *** '       $                              '*** MINIMIZATION FAILURE *** '
2875       $                              //'(mini2 in clouds_to_ctrack)'       $                              //'(clouds_to_ctrack)'
2876                                        print*,'initial guess: '
2877    
2878                                        print*,'AL_INI(1) = ',AL_INI(1)
2879                                        print*,'AL_INI(2) = ',AL_INI(2)
2880                                        print*,'AL_INI(3) = ',AL_INI(3)
2881                                        print*,'AL_INI(4) = ',AL_INI(4)
2882                                        print*,'AL_INI(5) = ',AL_INI(5)
2883                                   endif                                   endif
2884                                   chi2=-chi2  c                                 chi2=-chi2
2885                                endif                                endif
2886  *     **********************************************************  *     **********************************************************
2887  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
# Line 3035  c     $                                 Line 2894  c     $                                
2894  *     --------------------------  *     --------------------------
2895                                if(ntracks.eq.NTRACKSMAX)then                                if(ntracks.eq.NTRACKSMAX)then
2896                                                                    
2897                                   if(verbose)print*,                                   if(verbose.eq.1)print*,
2898       $                 '** warning ** number of candidate tracks '//       $                 '** warning ** number of candidate tracks '//
2899       $                 ' exceeds vector dimension '       $                 ' exceeds vector dimension '
2900       $                ,'( ',NTRACKSMAX,' )'       $                ,'( ',NTRACKSMAX,' )'
2901  c                                 good2=.false.  c                                 good2=.false.
2902  c                                 goto 880 !fill ntp and go to next event                      c                                 goto 880 !fill ntp and go to next event                    
2903                                     do iv=1,nviews
2904    c                                    mask_view(iv) = 7
2905                                        mask_view(iv) = mask_view(iv) + 2**6
2906                                     enddo
2907                                   iflag=1                                   iflag=1
2908                                   return                                   return
2909                                endif                                endif
2910                                                                
2911                                ntracks = ntracks + 1                                ntracks = ntracks + 1
2912                                                                
2913  c$$$                              ndof=0                                                                do ip=1,nplanes !top to bottom
2914                                do ip=1,nplanes  
 c$$$                                 ndof=ndof  
 c$$$     $                                +int(xgood(ip))  
 c$$$     $                                +int(ygood(ip))  
2915                                   XV_STORE(ip,ntracks)=sngl(xv(ip))                                   XV_STORE(ip,ntracks)=sngl(xv(ip))
2916                                   YV_STORE(ip,ntracks)=sngl(yv(ip))                                   YV_STORE(ip,ntracks)=sngl(yv(ip))
2917                                   ZV_STORE(ip,ntracks)=sngl(zv(ip))                                                                       ZV_STORE(ip,ntracks)=sngl(zv(ip))                                    
# Line 3067  c$$$     $                               Line 2927  c$$$     $                              
2927                                   AYV_STORE(ip,ntracks)=sngl(ayv(ip))                                   AYV_STORE(ip,ntracks)=sngl(ayv(ip))
2928                                   XGOOD_STORE(ip,ntracks)=sngl(xgood(ip))                                   XGOOD_STORE(ip,ntracks)=sngl(xgood(ip))
2929                                   YGOOD_STORE(ip,ntracks)=sngl(ygood(ip))                                   YGOOD_STORE(ip,ntracks)=sngl(ygood(ip))
2930    *                                NB! hit_plane is defined from bottom to top
2931                                   if(hit_plane(ip).ne.0)then                                   if(hit_plane(ip).ne.0)then
2932                                      CP_STORE(nplanes-ip+1,ntracks)=                                      CP_STORE(nplanes-ip+1,ntracks)=
2933       $                                   cp_match(ip,hit_plane(ip))       $                                   cp_match(ip,hit_plane(ip))
2934                                        SENSOR_STORE(nplanes-ip+1,ntracks)
2935         $                              = is_cp(cp_match(ip,hit_plane(ip)))
2936                                        LADDER_STORE(nplanes-ip+1,ntracks)
2937         $                                   = LADDER(
2938         $                                   clx(ip,icp_cp(
2939         $                                   cp_match(ip,hit_plane(ip)
2940         $                                   ))));
2941                                   else                                   else
2942                                      CP_STORE(nplanes-ip+1,ntracks)=0                                      CP_STORE(nplanes-ip+1,ntracks)=0
2943                                        SENSOR_STORE(nplanes-ip+1,ntracks)=0
2944                                        LADDER_STORE(nplanes-ip+1,ntracks)=0
2945                                   endif                                   endif
2946                                   CLS_STORE(nplanes-ip+1,ntracks)=0                                   BX_STORE(ip,ntracks)=0!I dont need it now
2947                                     BY_STORE(ip,ntracks)=0!I dont need it now
2948                                     CLS_STORE(ip,ntracks)=0
2949                                   do i=1,5                                   do i=1,5
2950                                      AL_STORE(i,ntracks)=sngl(AL(i))                                      AL_STORE(i,ntracks)=sngl(AL(i))
2951                                   enddo                                   enddo
2952                                enddo                                enddo
2953                                                                
 c$$$  *                             Number of Degree Of Freedom  
 c$$$  ndof=ndof-5                            
 c$$$  *                             reduced chi^2  
 c$$$  rchi2=chi2/dble(ndof)  
2954                                RCHI2_STORE(ntracks)=chi2                                RCHI2_STORE(ntracks)=chi2
2955                                                                
2956  *     --------------------------------  *     --------------------------------
# Line 3106  c$$$  rchi2=chi2/dble(ndof) Line 2974  c$$$  rchi2=chi2/dble(ndof)
2974           return           return
2975        endif        endif
2976                
2977        if(DEBUG)then  c$$$      if(DEBUG.EQ.1)then
2978           print*,'****** TRACK CANDIDATES ***********'  c$$$         print*,'****** TRACK CANDIDATES ***********'
2979           print*,'#         R. chi2        RIG'  c$$$         print*,'#         R. chi2        RIG'
2980           do i=1,ntracks  c$$$         do i=1,ntracks
2981              print*,i,' --- ',rchi2_store(i),' --- '  c$$$            print*,i,' --- ',rchi2_store(i),' --- '
2982       $           ,1./abs(AL_STORE(5,i))  c$$$     $           ,1./abs(AL_STORE(5,i))
2983           enddo  c$$$         enddo
2984           print*,'***********************************'  c$$$         print*,'***********************************'
2985    c$$$      endif
2986          if(DEBUG.EQ.1)then
2987            print*,'****** TRACK CANDIDATES *****************'
2988            print*,'#         R. chi2        RIG         ndof'
2989            do i=1,ntracks
2990              ndof=0                !(1)
2991              do ii=1,nplanes       !(1)
2992                ndof=ndof           !(1)
2993         $           +int(xgood_store(ii,i)) !(1)
2994         $           +int(ygood_store(ii,i)) !(1)
2995              enddo                 !(1)
2996              print*,i,' --- ',rchi2_store(i),' --- '
2997         $         ,1./abs(AL_STORE(5,i)),' --- ',ndof
2998            enddo
2999            print*,'*****************************************'
3000        endif        endif
3001                
3002                
# Line 3132  c$$$  rchi2=chi2/dble(ndof) Line 3015  c$$$  rchi2=chi2/dble(ndof)
3015    
3016        subroutine refine_track(ibest)        subroutine refine_track(ibest)
3017    
 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******************************************************  
3018    
3019        include 'commontracker.f'        include 'commontracker.f'
3020          include 'level1.f'
3021        include 'common_momanhough.f'        include 'common_momanhough.f'
3022        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3023        include 'common_mini_2.f'        include 'common_mini_2.f'
3024        include 'common_mech.f'        include 'common_mech.f'
       include 'momanhough_init.f'  
       include 'level1.f'  
3025        include 'calib.f'        include 'calib.f'
3026    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
3027  *     flag to chose PFA  *     flag to chose PFA
3028        character*10 PFA        character*10 PFA
3029        common/FINALPFA/PFA        common/FINALPFA/PFA
3030    
3031          real k(6)
3032          DATA k/1.099730,0.418900,0.220939,0.220907,0.418771,1.100674/
3033    
3034          real xp,yp,zp
3035          real xyzp(3),bxyz(3)
3036          equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3))
3037    
3038          if(DEBUG.EQ.1)print*,'refine_track:'
3039  *     =================================================  *     =================================================
3040  *     new estimate of positions using ETA algorithm  *     new estimate of positions using ETA algorithm
3041  *                          and  *                          and
# Line 3162  c      common/dbg/DEBUG Line 3044  c      common/dbg/DEBUG
3044        call track_init        call track_init
3045        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3046    
3047             xP=XV_STORE(nplanes-ip+1,ibest)
3048             yP=YV_STORE(nplanes-ip+1,ibest)
3049             zP=ZV_STORE(nplanes-ip+1,ibest)
3050             call gufld(xyzp,bxyz)
3051             BX_STORE(nplanes-ip+1,ibest)=bxyz(1)
3052             BY_STORE(nplanes-ip+1,ibest)=bxyz(2)
3053    c$$$  bxyz(1)=0
3054    c$$$         bxyz(2)=0
3055    c$$$         bxyz(3)=0
3056    *     |||||||||||||||||||||||||||||||||||||||||||||||||
3057  *     -------------------------------------------------  *     -------------------------------------------------
3058  *     If the plane has been already included, it just  *     If the plane has been already included, it just
3059  *     computes again the coordinates of the x-y couple  *     computes again the coordinates of the x-y couple
3060  *     using improved PFAs  *     using improved PFAs
3061  *     -------------------------------------------------  *     -------------------------------------------------
3062    *     |||||||||||||||||||||||||||||||||||||||||||||||||
3063           if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.           if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.
3064       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
3065                            
# Line 3180  c      common/dbg/DEBUG Line 3073  c      common/dbg/DEBUG
3073       $           ,ip_cp(id),ip       $           ,ip_cp(id),ip
3074              icx=clx(ip,icp)              icx=clx(ip,icp)
3075              icy=cly(ip,icp)              icy=cly(ip,icp)
3076    c            call xyz_PAM(icx,icy,is,
3077    c     $           PFA,PFA,
3078    c     $           AXV_STORE(nplanes-ip+1,ibest),
3079    c     $           AYV_STORE(nplanes-ip+1,ibest))
3080              call xyz_PAM(icx,icy,is,              call xyz_PAM(icx,icy,is,
 c     $           'ETA2','ETA2',  
3081       $           PFA,PFA,       $           PFA,PFA,
3082       $           AXV_STORE(nplanes-ip+1,ibest),       $           AXV_STORE(nplanes-ip+1,ibest),
3083       $           AYV_STORE(nplanes-ip+1,ibest))       $           AYV_STORE(nplanes-ip+1,ibest),
3084  c$$$  call xyz_PAM(icx,icy,is,       $           bxyz(1),
3085  c$$$  $              'COG2','COG2',       $           bxyz(2)
3086  c$$$  $              0.,       $           )
3087  c$$$  $              0.)  
3088              xm(nplanes-ip+1) = xPAM              xm(nplanes-ip+1) = xPAM
3089              ym(nplanes-ip+1) = yPAM              ym(nplanes-ip+1) = yPAM
3090              zm(nplanes-ip+1) = zPAM              zm(nplanes-ip+1) = zPAM
# Line 3197  c$$$  $              0.) Line 3093  c$$$  $              0.)
3093              resx(nplanes-ip+1) = resxPAM              resx(nplanes-ip+1) = resxPAM
3094              resy(nplanes-ip+1) = resyPAM              resy(nplanes-ip+1) = resyPAM
3095    
3096  c            dedxtrk(nplanes-ip+1) = (dedx(icx)+dedx(icy))/2. !(1)              dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx))
3097              dedxtrk_x(nplanes-ip+1)=dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)              dedxtrk_y(nplanes-ip+1)=sgnl(icy)/mip(VIEW(icy),LADDER(icy))
             dedxtrk_y(nplanes-ip+1)=dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)  
3098                            
3099    *     |||||||||||||||||||||||||||||||||||||||||||||||||
3100  *     -------------------------------------------------  *     -------------------------------------------------
3101  *     If the plane has NOT  been already included,  *     If the plane has NOT  been already included,
3102  *     it tries to include a COUPLE or a single cluster  *     it tries to include a COUPLE or a single cluster
3103  *     -------------------------------------------------  *     -------------------------------------------------
3104    *     |||||||||||||||||||||||||||||||||||||||||||||||||
3105           else                             else                  
3106                                
3107              xgood(nplanes-ip+1)=0              xgood(nplanes-ip+1)=0
# Line 3212  c            dedxtrk(nplanes-ip+1) = (de Line 3109  c            dedxtrk(nplanes-ip+1) = (de
3109                                
3110  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3111  *     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)  
3112              call whichsensor(ip,xP,yP,nldt,ist)              call whichsensor(ip,xP,yP,nldt,ist)
3113  *     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
3114              if(nldt.eq.0.or.ist.eq.0)goto 133              if(nldt.eq.0.or.ist.eq.0)goto 133
3115    
3116                SENSOR_STORE(nplanes-ip+1,IBEST)=ist
3117                LADDER_STORE(nplanes-ip+1,IBEST)=nldt
3118  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3119    
3120              if(DEBUG)then              if(DEBUG.EQ.1)then
3121                 print*,                 print*,
3122       $              '------ Plane ',ip,' intersected on LADDER ',nldt       $              '------ Plane ',ip,' intersected on LADDER ',nldt
3123       $              ,' SENSOR ',ist       $              ,' SENSOR ',ist
# Line 3231  c            dedxtrk(nplanes-ip+1) = (de Line 3128  c            dedxtrk(nplanes-ip+1) = (de
3128  *     ===========================================  *     ===========================================
3129  *     STEP 1 >>>>>>>  try to include a new couple  *     STEP 1 >>>>>>>  try to include a new couple
3130  *     ===========================================  *     ===========================================
3131  c            if(DEBUG)print*,'>>>> try to include a new couple'  c            if(DEBUG.EQ.1)print*,'>>>> try to include a new couple'
3132              distmin=1000000.              distmin=1000000.
3133              xmm = 0.              xmm = 0.
3134              ymm = 0.              ymm = 0.
# Line 3254  c     $              cl_used(icy).eq.1.o Line 3151  c     $              cl_used(icy).eq.1.o
3151  *            *          
3152                 call xyz_PAM(icx,icy,ist,                 call xyz_PAM(icx,icy,ist,
3153       $              PFA,PFA,       $              PFA,PFA,
 c     $              'ETA2','ETA2',  
3154       $              AXV_STORE(nplanes-ip+1,ibest),       $              AXV_STORE(nplanes-ip+1,ibest),
3155       $              AYV_STORE(nplanes-ip+1,ibest))       $              AYV_STORE(nplanes-ip+1,ibest),
3156         $              bxyz(1),
3157         $              bxyz(2)
3158         $              )
3159                                
3160                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3161    c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3162                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
3163                 if(DEBUG)print*,'( couple ',id                 if(DEBUG.EQ.1)print*,'( couple ',id
3164       $              ,' ) normalized distance ',distance       $              ,' ) distance ',distance
3165                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3166                    xmm = xPAM                    xmm = xPAM
3167                    ymm = yPAM                    ymm = yPAM
# Line 3270  c     $              'ETA2','ETA2', Line 3170  c     $              'ETA2','ETA2',
3170                    rymm = resyPAM                    rymm = resyPAM
3171                    distmin = distance                    distmin = distance
3172                    idm = id                                      idm = id                  
3173  c                 dedxmm = (dedx(icx)+dedx(icy))/2. !(1)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3174                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3175                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)  c     QUIQUI --> non devo moltiplicare per clinc?!?!?!
3176                      clincnewc=10*sqrt(rymm**2+rxmm**2 !QUIQUI
3177         $                 +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2))) !QUIQUI
3178                 endif                 endif
3179   1188          continue   1188          continue
3180              enddo               !end loop on couples on plane icp              enddo               !end loop on couples on plane icp
3181              if(distmin.le.clinc)then                    c            if(distmin.le.clinc)then     !QUIQUI              
3182                if(distmin.le.clincnewc)then     !QUIQUI              
3183  *              -----------------------------------  *              -----------------------------------
3184                 xm(nplanes-ip+1) = xmm         !<<<                 xm(nplanes-ip+1) = xmm !<<<
3185                 ym(nplanes-ip+1) = ymm         !<<<                 ym(nplanes-ip+1) = ymm !<<<
3186                 zm(nplanes-ip+1) = zmm         !<<<                 zm(nplanes-ip+1) = zmm !<<<
3187                 xgood(nplanes-ip+1) = 1        !<<<                 xgood(nplanes-ip+1) = 1 !<<<
3188                 ygood(nplanes-ip+1) = 1        !<<<                 ygood(nplanes-ip+1) = 1 !<<<
3189                 resx(nplanes-ip+1)=rxmm        !<<<                 resx(nplanes-ip+1)=rxmm !<<<
3190                 resy(nplanes-ip+1)=rymm        !<<<                 resy(nplanes-ip+1)=rymm !<<<
3191  c              dedxtrk(nplanes-ip+1) = dedxmm !<<<  !(1)                 dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3192                 dedxtrk_x(nplanes-ip+1) = dedxmmx    !(1)                 dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
                dedxtrk_y(nplanes-ip+1) = dedxmmy    !(1)  
3193  *              -----------------------------------  *              -----------------------------------
3194                 CP_STORE(nplanes-ip+1,ibest)=idm                       CP_STORE(nplanes-ip+1,ibest)=idm      
3195                 if(DEBUG)print*,'%%%% included couple ',idm                 if(DEBUG.EQ.1)print*,'%%%% included couple ',idm
3196       $              ,' (norm.dist.= ',distmin,', cut ',clinc,' )'       $              ,' (dist.= ',distmin,', cut ',clinc,' )'
3197                 goto 133         !next plane                 goto 133         !next plane
3198              endif              endif
3199  *     ================================================  *     ================================================
3200  *     STEP 2 >>>>>>>  try to include a single cluster  *     STEP 2 >>>>>>>  try to include a single cluster
3201  *                     either from a couple or single  *                     either from a couple or single
3202  *     ================================================  *     ================================================
3203  c            if(DEBUG)print*,'>>>> try to include a new cluster'  c            if(DEBUG.EQ.1)print*,'>>>> try to include a new cluster'
3204              distmin=1000000.              distmin=1000000.
3205              xmm_A = 0.          !---------------------------              xmm_A = 0.          !---------------------------
3206              ymm_A = 0.          ! init variables that              ymm_A = 0.          ! init variables that
# Line 3324  c            if(DEBUG)print*,'>>>> try t Line 3226  c            if(DEBUG)print*,'>>>> try t
3226  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
3227                 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)
3228  *                                              !jump to the Y cluster  *                                              !jump to the Y cluster
3229    c               call xyz_PAM(icx,0,ist,
3230    c     $              PFA,PFA,
3231    c     $              AXV_STORE(nplanes-ip+1,ibest),0.)              
3232                 call xyz_PAM(icx,0,ist,                 call xyz_PAM(icx,0,ist,
 c     $              'ETA2','ETA2',  
3233       $              PFA,PFA,       $              PFA,PFA,
3234       $              AXV_STORE(nplanes-ip+1,ibest),0.)                     $              AXV_STORE(nplanes-ip+1,ibest),0.,
3235         $              bxyz(1),
3236         $              bxyz(2)
3237         $              )              
3238                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3239  c     if(DEBUG)print*,'normalized distance ',distance  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3240                 if(DEBUG)print*,'( cl-X ',icx                 if(DEBUG.EQ.1)print*,'( cl-X ',icx
3241       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3242                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3243                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3244                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 3343  c     if(DEBUG)print*,'normalized distan Line 3250  c     if(DEBUG)print*,'normalized distan
3250                    rymm = resyPAM                    rymm = resyPAM
3251                    distmin = distance                    distmin = distance
3252                    iclm = icx                    iclm = icx
3253  c                  dedxmm = dedx(icx) !(1)  c                  dedxmm = sgnl(icx) !(1)
3254                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3255                    dedxmmy = 0.        !(1)                    dedxmmy = 0.        !(1)
3256                 endif                                   endif                  
3257  11881          continue  11881          continue
# Line 3352  c                  dedxmm = dedx(icx) !( Line 3259  c                  dedxmm = dedx(icx) !(
3259  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
3260                 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)
3261  *                                              !jump to the next couple  *                                              !jump to the next couple
3262    c               call xyz_PAM(0,icy,ist,
3263    c     $              PFA,PFA,
3264    c     $              0.,AYV_STORE(nplanes-ip+1,ibest))
3265                 call xyz_PAM(0,icy,ist,                 call xyz_PAM(0,icy,ist,
 c     $              'ETA2','ETA2',  
3266       $              PFA,PFA,       $              PFA,PFA,
3267       $              0.,AYV_STORE(nplanes-ip+1,ibest))       $              0.,AYV_STORE(nplanes-ip+1,ibest),
3268         $              bxyz(1),
3269         $              bxyz(2)
3270         $              )
3271                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3272                 if(DEBUG)print*,'( cl-Y ',icy  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3273       $              ,' in cp ',id,' ) normalized distance ',distance                 if(DEBUG.EQ.1)print*,'( cl-Y ',icy
3274         $              ,' in cp ',id,' ) distance ',distance
3275                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3276                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3277                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 3370  c     $              'ETA2','ETA2', Line 3283  c     $              'ETA2','ETA2',
3283                    rymm = resyPAM                    rymm = resyPAM
3284                    distmin = distance                    distmin = distance
3285                    iclm = icy                    iclm = icy
3286  c                 dedxmm = dedx(icy)  !(1)  c                 dedxmm = sgnl(icy)  !(1)
3287                    dedxmmx = 0.        !(1)                    dedxmmx = 0.        !(1)
3288                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3289                 endif                                   endif                  
3290  11882          continue  11882          continue
3291              enddo               !end loop on cluster inside couples              enddo               !end loop on cluster inside couples
3292  *----- single clusters -----------------------------------------------      *----- single clusters -----------------------------------------------  
3293    c            print*,'## ncls(',ip,') ',ncls(ip)
3294              do ic=1,ncls(ip)    !loop on single clusters              do ic=1,ncls(ip)    !loop on single clusters
3295                 icl=cls(ip,ic)                 icl=cls(ip,ic)
3296  c               if(cl_used(icl).eq.1.or.     !if the cluster is already used  c               if(cl_used(icl).eq.1.or.     !if the cluster is already used
# Line 3385  c               if(cl_used(icl).eq.1.or. Line 3299  c               if(cl_used(icl).eq.1.or.
3299       $              .false.)goto 18882      !jump to the next singlet       $              .false.)goto 18882      !jump to the next singlet
3300                 if(mod(VIEW(icl),2).eq.0)then!<---- X view                 if(mod(VIEW(icl),2).eq.0)then!<---- X view
3301                    call xyz_PAM(icl,0,ist,                    call xyz_PAM(icl,0,ist,
 c     $                 'ETA2','ETA2',  
3302       $                 PFA,PFA,       $                 PFA,PFA,
3303       $                 AXV_STORE(nplanes-ip+1,ibest),0.)       $                 AXV_STORE(nplanes-ip+1,ibest),0.,
3304         $                 bxyz(1),
3305         $                 bxyz(2)
3306         $                 )
3307                 else                         !<---- Y view                 else                         !<---- Y view
3308                    call xyz_PAM(0,icl,ist,                    call xyz_PAM(0,icl,ist,
 c     $                 'ETA2','ETA2',  
3309       $                 PFA,PFA,       $                 PFA,PFA,
3310       $                 0.,AYV_STORE(nplanes-ip+1,ibest))       $                 0.,AYV_STORE(nplanes-ip+1,ibest),
3311         $                 bxyz(1),
3312         $                 bxyz(2)
3313         $                 )
3314                 endif                 endif
3315    
3316                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3317                 if(DEBUG)print*,'( cl-s ',icl  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3318       $              ,' ) normalized distance ',distance                 if(DEBUG.EQ.1)print*,'( cl-s ',icl
3319         $              ,' ) distance ',distance
3320                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3321    c                  if(DEBUG.EQ.1)print*,'YES'
3322                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3323                    ymm_A = yPAM_A                    ymm_A = yPAM_A
3324                    zmm_A = zPAM_A                    zmm_A = zPAM_A
# Line 3409  c     $                 'ETA2','ETA2', Line 3329  c     $                 'ETA2','ETA2',
3329                    rymm = resyPAM                    rymm = resyPAM
3330                    distmin = distance                      distmin = distance  
3331                    iclm = icl                    iclm = icl
 c                  dedxmm = dedx(icl)                   !(1)  
3332                    if(mod(VIEW(icl),2).eq.0)then !<---- X view                    if(mod(VIEW(icl),2).eq.0)then !<---- X view
3333                       dedxmmx = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmx = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3334                       dedxmmy = 0.                       !(1)                       dedxmmy = 0.                  
3335                    else          !<---- Y view                    else          !<---- Y view
3336                       dedxmmx = 0.                       !(1)                       dedxmmx = 0.                  
3337                       dedxmmy = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmy = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3338                    endif                    endif
3339                 endif                                   endif                  
3340  18882          continue  18882          continue
3341              enddo               !end loop on single clusters              enddo               !end loop on single clusters
3342    c            print*,'## distmin ', distmin,' clinc ',clinc
3343    
3344              if(distmin.le.clinc)then                    c     QUIQUI------------
3345                  c     anche qui: non ci vuole clinc???
3346                 CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<                  if(iclm.ne.0)then
 *              ----------------------------  
3347                 if(mod(VIEW(iclm),2).eq.0)then                 if(mod(VIEW(iclm),2).eq.0)then
3348                    XGOOD(nplanes-ip+1)=1.                    clincnew=
3349                    resx(nplanes-ip+1)=rxmm       $                 20*
3350                    if(DEBUG)print*,'%%%% included X-cl ',iclm       $                 sqrt(rxmm**2+RCHI2_STORE(ibest)*k(ip)*cov(1,1))
3351       $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'                 else if(mod(VIEW(iclm),2).ne.0)then
3352                 else                    clincnew=
3353                    YGOOD(nplanes-ip+1)=1.       $                 10*
3354                    resy(nplanes-ip+1)=rymm       $                 sqrt(rymm**2+RCHI2_STORE(ibest)*k(ip)*cov(2,2))
3355                    if(DEBUG)print*,'%%%% included Y-cl ',iclm                 endif
3356       $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'  c     QUIQUI------------
3357                  
3358                   if(distmin.le.clincnew)then   !QUIQUI
3359    c     if(distmin.le.clinc)then          !QUIQUI          
3360                      
3361                      CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<    
3362    *     ----------------------------
3363    c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
3364                      if(mod(VIEW(iclm),2).eq.0)then
3365                         XGOOD(nplanes-ip+1)=1.
3366                         resx(nplanes-ip+1)=rxmm
3367                         if(DEBUG.EQ.1)print*,'%%%% included X-cl ',iclm
3368         $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3369         $                    ,', dist.= ',distmin
3370         $                    ,', cut ',clinc,' )'
3371                      else
3372                         YGOOD(nplanes-ip+1)=1.
3373                         resy(nplanes-ip+1)=rymm
3374                         if(DEBUG.EQ.1)print*,'%%%% included Y-cl ',iclm
3375         $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3376         $                    ,', dist.= ', distmin
3377         $                    ,', cut ',clinc,' )'
3378                      endif
3379    c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
3380    *     ----------------------------
3381                      xm_A(nplanes-ip+1) = xmm_A
3382                      ym_A(nplanes-ip+1) = ymm_A
3383                      xm_B(nplanes-ip+1) = xmm_B
3384                      ym_B(nplanes-ip+1) = ymm_B
3385                      zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.
3386                      dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3387                      dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
3388    *     ----------------------------
3389                 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)  
 *              ----------------------------  
3390              endif              endif
3391           endif           endif
3392   133     continue   133     continue
# Line 3464  c              dedxtrk(nplanes-ip+1) = d Line 3405  c              dedxtrk(nplanes-ip+1) = d
3405  *                                                 *  *                                                 *
3406  *                                                 *  *                                                 *
3407  **************************************************  **************************************************
 cccccc 12/08/2006 modified by elena ---> (1)  
3408  *  *
3409        subroutine clean_XYclouds(ibest,iflag)        subroutine clean_XYclouds(ibest,iflag)
3410    
3411        include 'commontracker.f'        include 'commontracker.f'
3412          include 'level1.f'
3413        include 'common_momanhough.f'        include 'common_momanhough.f'
3414        include 'momanhough_init.f'        include 'level2.f'      
       include 'level2.f'        !(1)  
 c      include 'calib.f'  
 c      include 'level1.f'  
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
3415    
3416          if(DEBUG.EQ.1)print*,'clean_XYclouds:'
3417    
3418        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3419    
# Line 3487  c      common/dbg/DEBUG Line 3423  c      common/dbg/DEBUG
3423              if(id.ne.0)then              if(id.ne.0)then
3424                 iclx=clx(ip,icp_cp(id))                 iclx=clx(ip,icp_cp(id))
3425                 icly=cly(ip,icp_cp(id))                 icly=cly(ip,icp_cp(id))
3426  c               cl_used(iclx)=1  !tag used clusters  c$$$               cl_used(iclx)=ntrk  !tag used clusters
3427  c               cl_used(icly)=1  !tag used clusters  c$$$               cl_used(icly)=ntrk  !tag used clusters
                cl_used(iclx)=ntrk  !tag used clusters !(1)  
                cl_used(icly)=ntrk  !tag used clusters !(1)  
3428              elseif(icl.ne.0)then              elseif(icl.ne.0)then
3429  c               cl_used(icl)=1   !tag used clusters  c$$$               cl_used(icl)=ntrk   !tag used clusters
                cl_used(icl)=ntrk   !tag used clusters !1)  
3430              endif              endif
3431                            
 c               if(DEBUG)then  
 c                  print*,ip,' <<< ',id  
 c               endif  
3432  *     -----------------------------  *     -----------------------------
3433  *     remove the couple from clouds  *     remove the couple from clouds
3434  *     remove also vitual couples containing the  *     remove also vitual couples containing the
# Line 3515  c               endif Line 3445  c               endif
3445       $              cly(ip,icp).eq.icl       $              cly(ip,icp).eq.icl
3446       $              )then       $              )then
3447                    id=id_cp(ip,icp,1)                    id=id_cp(ip,icp,1)
3448                    if(DEBUG)then                    if(DEBUG.EQ.1)then
3449                       print*,ip,' <<< cp ',id                       print*,ip,' <<< cp ',id
3450       $                    ,' ( cl-x '       $                    ,' ( cl-x '
3451       $                    ,clx(ip,icp)       $                    ,clx(ip,icp)
# Line 3549  c               endif Line 3479  c               endif
3479    
3480    
3481    
 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$$$  
   
3482    
3483    
3484  *     ****************************************************  *     ****************************************************
3485    
3486        subroutine init_level2        subroutine init_level2
3487    
 c*****************************************************  
 c     07/10/2005 modified by elena vannuccini --> (1)  
 c*****************************************************  
   
3488        include 'commontracker.f'        include 'commontracker.f'
3489          include 'level1.f'
3490        include 'common_momanhough.f'        include 'common_momanhough.f'
3491        include 'level2.f'        include 'level2.f'
       include 'level1.f'  
3492    
3493    *     ---------------------------------
3494    *     variables initialized from level1
3495    *     ---------------------------------
3496        do i=1,nviews        do i=1,nviews
3497           good2(i)=good1(i)           good2(i)=good1(i)
3498             do j=1,nva1_view
3499                vkflag(i,j)=1
3500                if(cnnev(i,j).le.0)then
3501                   vkflag(i,j)=cnnev(i,j)
3502                endif
3503             enddo
3504        enddo        enddo
3505    *     ----------------
3506  c      good2 = 0!.false.  *     level2 variables
3507  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*****************************************************  
   
3508        NTRK = 0        NTRK = 0
3509        do it=1,NTRKMAX!NTRACKSMAX        do it=1,NTRKMAX
3510           IMAGE(IT)=0           IMAGE(IT)=0
3511           CHI2_nt(IT) = -100000.           CHI2_nt(IT) = -100000.
 c         BdL(IT) = 0.  
3512           do ip=1,nplanes           do ip=1,nplanes
3513              XM_nt(IP,IT) = 0              XM_nt(IP,IT) = 0
3514              YM_nt(IP,IT) = 0              YM_nt(IP,IT) = 0
3515              ZM_nt(IP,IT) = 0              ZM_nt(IP,IT) = 0
3516              RESX_nt(IP,IT) = 0              RESX_nt(IP,IT) = 0
3517              RESY_nt(IP,IT) = 0              RESY_nt(IP,IT) = 0
3518                TAILX_nt(IP,IT) = 0
3519                TAILY_nt(IP,IT) = 0
3520                XBAD(IP,IT) = 0
3521                YBAD(IP,IT) = 0
3522              XGOOD_nt(IP,IT) = 0              XGOOD_nt(IP,IT) = 0
3523              YGOOD_nt(IP,IT) = 0              YGOOD_nt(IP,IT) = 0
3524  c*****************************************************              LS(IP,IT) = 0
 cccccc 11/9/2005 modified by david fedele  
3525              DEDX_X(IP,IT) = 0              DEDX_X(IP,IT) = 0
3526              DEDX_Y(IP,IT) = 0              DEDX_Y(IP,IT) = 0
 c******************************************************  
 cccccc 17/8/2006 modified by elena  
3527              CLTRX(IP,IT) = 0              CLTRX(IP,IT) = 0
3528              CLTRY(IP,IT) = 0              CLTRY(IP,IT) = 0
3529                multmaxx(ip,it) = 0
3530                seedx(ip,it)    = 0
3531                xpu(ip,it)      = 0
3532                multmaxy(ip,it) = 0
3533                seedy(ip,it)    = 0
3534                ypu(ip,it)      = 0
3535           enddo           enddo
3536           do ipa=1,5           do ipa=1,5
3537              AL_nt(IPA,IT) = 0              AL_nt(IPA,IT) = 0
# Line 3717  cccccc 17/8/2006 modified by elena Line 3540  cccccc 17/8/2006 modified by elena
3540              enddo                                enddo                  
3541           enddo                             enddo                  
3542        enddo        enddo
         
         
 c*****************************************************  
 cccccc 11/9/2005 modified by david fedele  
3543        nclsx=0        nclsx=0
3544        nclsy=0              nclsy=0      
3545        do ip=1,NSINGMAX        do ip=1,NSINGMAX
3546          planex(ip)=0          planex(ip)=0
 c        xs(ip)=0  
3547          xs(1,ip)=0          xs(1,ip)=0
3548          xs(2,ip)=0          xs(2,ip)=0
3549          sgnlxs(ip)=0          sgnlxs(ip)=0
3550          planey(ip)=0          planey(ip)=0
 c        ys(ip)=0  
3551          ys(1,ip)=0          ys(1,ip)=0
3552          ys(2,ip)=0          ys(2,ip)=0
3553          sgnlys(ip)=0          sgnlys(ip)=0
3554        enddo        enddo
 c*******************************************************  
3555        end        end
3556    
3557    
# Line 3750  c*************************************** Line 3566  c***************************************
3566  ************************************************************  ************************************************************
3567    
3568    
3569          subroutine init_hough
3570    
3571          include 'commontracker.f'
3572          include 'level1.f'
3573          include 'common_momanhough.f'
3574          include 'common_hough.f'
3575          include 'level2.f'
3576    
3577          ntrpt_nt=0
3578          ndblt_nt=0
3579          NCLOUDS_XZ_nt=0
3580          NCLOUDS_YZ_nt=0
3581          do idb=1,ndblt_max_nt
3582             db_cloud_nt(idb)=0
3583             alfayz1_nt(idb)=0      
3584             alfayz2_nt(idb)=0      
3585          enddo
3586          do itr=1,ntrpt_max_nt
3587             tr_cloud_nt(itr)=0
3588             alfaxz1_nt(itr)=0      
3589             alfaxz2_nt(itr)=0      
3590             alfaxz3_nt(itr)=0      
3591          enddo
3592          do idb=1,ncloyz_max      
3593            ptcloud_yz_nt(idb)=0    
3594            alfayz1_av_nt(idb)=0    
3595            alfayz2_av_nt(idb)=0    
3596          enddo                    
3597          do itr=1,ncloxz_max      
3598            ptcloud_xz_nt(itr)=0    
3599            alfaxz1_av_nt(itr)=0    
3600            alfaxz2_av_nt(itr)=0    
3601            alfaxz3_av_nt(itr)=0    
3602          enddo                    
3603    
3604          ntrpt=0                  
3605          ndblt=0                  
3606          NCLOUDS_XZ=0              
3607          NCLOUDS_YZ=0              
3608          do idb=1,ndblt_max        
3609            db_cloud(idb)=0        
3610            cpyz1(idb)=0            
3611            cpyz2(idb)=0            
3612            alfayz1(idb)=0          
3613            alfayz2(idb)=0          
3614          enddo                    
3615          do itr=1,ntrpt_max        
3616            tr_cloud(itr)=0        
3617            cpxz1(itr)=0            
3618            cpxz2(itr)=0            
3619            cpxz3(itr)=0            
3620            alfaxz1(itr)=0          
3621            alfaxz2(itr)=0          
3622            alfaxz3(itr)=0          
3623          enddo                    
3624          do idb=1,ncloyz_max      
3625            ptcloud_yz(idb)=0      
3626            alfayz1_av(idb)=0      
3627            alfayz2_av(idb)=0      
3628            do idbb=1,ncouplemaxtot
3629              cpcloud_yz(idb,idbb)=0
3630            enddo                  
3631          enddo                    
3632          do itr=1,ncloxz_max      
3633            ptcloud_xz(itr)=0      
3634            alfaxz1_av(itr)=0      
3635            alfaxz2_av(itr)=0      
3636            alfaxz3_av(itr)=0      
3637            do itrr=1,ncouplemaxtot
3638              cpcloud_xz(itr,itrr)=0
3639            enddo                  
3640          enddo                    
3641          end
3642    ************************************************************
3643    *
3644    *
3645    *
3646    *
3647    *
3648    *
3649    *
3650    ************************************************************
3651    
3652    
3653        subroutine fill_level2_tracks(ntr)        subroutine fill_level2_tracks(ntr)
3654    
3655  *     -------------------------------------------------------  *     -------------------------------------------------------
# Line 3761  c*************************************** Line 3661  c***************************************
3661            
3662        include 'commontracker.f'        include 'commontracker.f'
3663        include 'level1.f'        include 'level1.f'
3664          include 'common_momanhough.f'
3665        include 'level2.f'        include 'level2.f'
3666        include 'common_mini_2.f'        include 'common_mini_2.f'
3667        include 'common_momanhough.f'        include 'calib.f'
3668        real sinth,phi,pig        !(4)  
3669          character*10 PFA
3670          common/FINALPFA/PFA
3671    
3672          real sinth,phi,pig
3673          integer ssensor,sladder
3674        pig=acos(-1.)        pig=acos(-1.)
3675    
3676  c      good2=1!.true.  *     -------------------------------------
3677        chi2_nt(ntr)        = sngl(chi2)        chi2_nt(ntr)        = sngl(chi2)
3678        nstep_nt(ntr)       = nstep        nstep_nt(ntr)       = nstep
3679    *     -------------------------------------
3680        phi   = al(4)             !(4)        phi   = al(4)          
3681        sinth = al(3)             !(4)        sinth = al(3)            
3682        if(sinth.lt.0)then        !(4)        if(sinth.lt.0)then      
3683           sinth = -sinth         !(4)           sinth = -sinth        
3684           phi = phi + pig        !(4)           phi = phi + pig      
3685        endif                     !(4)        endif                    
3686        npig = aint(phi/(2*pig))  !(4)        npig = aint(phi/(2*pig))
3687        phi = phi - npig*2*pig    !(4)        phi = phi - npig*2*pig  
3688        if(phi.lt.0)              !(4)        if(phi.lt.0)            
3689       $     phi = phi + 2*pig    !(4)       $     phi = phi + 2*pig  
3690        al(4) = phi               !(4)        al(4) = phi              
3691        al(3) = sinth             !(4)        al(3) = sinth            
 *****************************************************  
3692        do i=1,5        do i=1,5
3693           al_nt(i,ntr)     = sngl(al(i))           al_nt(i,ntr)     = sngl(al(i))
3694           do j=1,5           do j=1,5
3695              coval(i,j,ntr) = sngl(cov(i,j))              coval(i,j,ntr) = sngl(cov(i,j))
3696           enddo           enddo
 c     print*,al_nt(i,ntr)  
3697        enddo        enddo
3698          *     -------------------------------------      
3699        do ip=1,nplanes           ! loop on planes        do ip=1,nplanes           ! loop on planes
3700           xgood_nt(ip,ntr) = int(xgood(ip))           xgood_nt(ip,ntr) = int(xgood(ip))
3701           ygood_nt(ip,ntr) = int(ygood(ip))           ygood_nt(ip,ntr) = int(ygood(ip))
# Line 3800  c     print*,al_nt(i,ntr) Line 3704  c     print*,al_nt(i,ntr)
3704           zm_nt(ip,ntr)    = sngl(zm(ip))           zm_nt(ip,ntr)    = sngl(zm(ip))
3705           RESX_nt(IP,ntr)  = sngl(resx(ip))           RESX_nt(IP,ntr)  = sngl(resx(ip))
3706           RESY_nt(IP,ntr)  = sngl(resy(ip))           RESY_nt(IP,ntr)  = sngl(resy(ip))
3707             TAILX_nt(IP,ntr) = 0.
3708             TAILY_nt(IP,ntr) = 0.
3709           xv_nt(ip,ntr)    = sngl(xv(ip))           xv_nt(ip,ntr)    = sngl(xv(ip))
3710           yv_nt(ip,ntr)    = sngl(yv(ip))           yv_nt(ip,ntr)    = sngl(yv(ip))
3711           zv_nt(ip,ntr)    = sngl(zv(ip))           zv_nt(ip,ntr)    = sngl(zv(ip))
3712           axv_nt(ip,ntr)   = sngl(axv(ip))           axv_nt(ip,ntr)   = sngl(axv(ip))
3713           ayv_nt(ip,ntr)   = sngl(ayv(ip))           ayv_nt(ip,ntr)   = sngl(ayv(ip))  
3714  c        dedxp(ip,ntr)    = sngl(dedxtrk(ip))   !(1)  
3715           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)           factor = sqrt(
3716           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)         $        tan( acos(-1.) * sngl(axv(ip)) /180. )**2 +
3717         $        tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 +
3718         $        1. )
3719    
3720             dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)
3721             dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)  
3722        
3723           id  = CP_STORE(ip,IDCAND)           ax   = axv_nt(ip,ntr)
3724             ay   = ayv_nt(ip,ntr)
3725             bfx  = BX_STORE(ip,IDCAND)
3726             bfy  = BY_STORE(ip,IDCAND)
3727    c$$$         if(ip.eq.6) ax = -1. * axv_nt(ip,ntr)
3728    c$$$         if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND)
3729    c$$$         tgtemp   = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001
3730    c$$$         angx     = 180.*atan(tgtemp)/acos(-1.)
3731    c$$$         tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
3732    c$$$         angy    = 180.*atan(tgtemp)/acos(-1.)
3733    
3734             angx = effectiveangle(ax,2*ip,bfy)
3735             angy = effectiveangle(ay,2*ip-1,bfx)
3736            
3737            
3738    c         print*,'* ',ip,bfx,bfy,angx,angy
3739    
3740             id  = CP_STORE(ip,IDCAND) ! couple id
3741           icl = CLS_STORE(ip,IDCAND)           icl = CLS_STORE(ip,IDCAND)
3742             ssensor = -1
3743             sladder = -1
3744             ssensor = SENSOR_STORE(ip,IDCAND)
3745             sladder = LADDER_STORE(ip,IDCAND)
3746             if(ip.eq.6.and.ssensor.ne.0)ssensor = 3 - ssensor !notazione paolo x align
3747             LS(IP,ntr)      = ssensor+10*sladder
3748    
3749           if(id.ne.0)then           if(id.ne.0)then
3750    c           >>> is a couple
3751              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))
3752              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))
3753  c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)  
3754                cl_used(cltrx(ip,ntr)) = 1     !tag used clusters          
3755                cl_used(cltry(ip,ntr)) = 1     !tag used clusters          
3756    
3757                xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))
3758                ybad(ip,ntr)= nbadstrips(4,cly(nplanes-ip+1,icp_cp(id)))
3759    
3760    
3761                if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0)
3762         $           dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3763                if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)
3764         $           dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3765    
3766                multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
3767         $                         +10000*mult(cltrx(ip,ntr))
3768                seedx(ip,ntr)    = clsignal(indmax(cltrx(ip,ntr)))
3769         $           /clsigma(indmax(cltrx(ip,ntr)))
3770                call applypfa(PFA,cltrx(ip,ntr),angx,corr,res)
3771                xpu(ip,ntr)      = corr
3772    
3773                multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
3774         $                         +10000*mult(cltry(ip,ntr))
3775                seedy(ip,ntr)    = clsignal(indmax(cltry(ip,ntr)))
3776         $           /clsigma(indmax(cltry(ip,ntr)))
3777                call applypfa(PFA,cltry(ip,ntr),angy,corr,res)
3778                ypu(ip,ntr)      = corr
3779    
3780           elseif(icl.ne.0)then           elseif(icl.ne.0)then
3781              if(mod(VIEW(icl),2).eq.0)cltrx(ip,ntr)=icl  
3782              if(mod(VIEW(icl),2).eq.1)cltry(ip,ntr)=icl              cl_used(icl) = 1    !tag used clusters          
3783  c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)  
3784                if(mod(VIEW(icl),2).eq.0)then
3785                   cltrx(ip,ntr)=icl
3786    
3787                   xbad(ip,ntr) = nbadstrips(4,icl)
3788    
3789                   if(nsatstrips(icl).gt.0)dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3790    
3791                   multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
3792         $                         +10000*mult(cltrx(ip,ntr))
3793                   seedx(ip,ntr)    = clsignal(indmax(cltrx(ip,ntr)))
3794         $           /clsigma(indmax(cltrx(ip,ntr)))
3795                   call applypfa(PFA,cltrx(ip,ntr),angx,corr,res)
3796                   xpu(ip,ntr)      = corr
3797    
3798                elseif(mod(VIEW(icl),2).eq.1)then
3799                   cltry(ip,ntr)=icl
3800    
3801                   ybad(ip,ntr) = nbadstrips(4,icl)
3802    
3803                   if(nsatstrips(icl).gt.0)dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3804    
3805                   multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
3806         $                         +10000*mult(cltry(ip,ntr))
3807                   seedy(ip,ntr)    = clsignal(indmax(cltry(ip,ntr)))
3808         $           /clsigma(indmax(cltry(ip,ntr)))
3809                   call applypfa(PFA,cltry(ip,ntr),angy,corr,res)
3810                   ypu(ip,ntr)      = corr
3811                  
3812                endif
3813    
3814           endif                     endif          
3815    
3816        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)  
3817    
3818          if(DEBUG.eq.1)then
3819             print*,'> STORING TRACK ',ntr
3820             print*,'clusters: '
3821             do ip=1,6
3822                print*,'> ',ip,' -- ',cltrx(ip,ntr),cltry(ip,ntr)
3823             enddo
3824          endif
3825    
3826    c$$$      print*,(xgood(i),i=1,6)
3827    c$$$      print*,(ygood(i),i=1,6)
3828    c$$$      print*,(ls(i,ntr),i=1,6)
3829    c$$$      print*,(dedx_x(i,ntr),i=1,6)
3830    c$$$      print*,(dedx_y(i,ntr),i=1,6)
3831    c$$$      print*,'-----------------------'
3832    
3833        end        end
3834    
3835        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*****************************************************  
3836    
3837  *     -------------------------------------------------------  *     -------------------------------------------------------
3838  *     This routine fills the  elements of the variables  *     This routine fills the  elements of the variables
# Line 3846  c*************************************** Line 3841  c***************************************
3841  *     -------------------------------------------------------  *     -------------------------------------------------------
3842    
3843        include 'commontracker.f'        include 'commontracker.f'
       include 'level1.f'  
       include 'level2.f'  
3844        include 'calib.f'        include 'calib.f'
3845          include 'level1.f'
3846        include 'common_momanhough.f'        include 'common_momanhough.f'
3847          include 'level2.f'
3848        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3849    
3850  *     count #cluster per plane not associated to any track  *     count #cluster per plane not associated to any track
 c      good2=1!.true.  
3851        nclsx = 0        nclsx = 0
3852        nclsy = 0        nclsy = 0
3853    
3854          do iv = 1,nviews
3855    c         if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv)
3856             good2(iv) = good2(iv) + mask_view(iv)*2**8
3857          enddo
3858    
3859          if(DEBUG.eq.1)then
3860             print*,'> STORING SINGLETS '
3861          endif
3862    
3863        do icl=1,nclstr1        do icl=1,nclstr1
3864    
3865             ip=nplanes-npl(VIEW(icl))+1            
3866            
3867           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
             ip=nplanes-npl(VIEW(icl))+1              
3868              if(mod(VIEW(icl),2).eq.0)then !=== X views              if(mod(VIEW(icl),2).eq.0)then !=== X views
3869                 nclsx = nclsx + 1                 nclsx = nclsx + 1
3870                 planex(nclsx) = ip                 planex(nclsx) = ip
3871                 sgnlxs(nclsx) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3872                   if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx)
3873                 clsx(nclsx)   = icl                 clsx(nclsx)   = icl
3874                 do is=1,2                 do is=1,2
3875  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
3876                    call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)  c                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
3877                      call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.,0.,0.)
3878                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2
3879                 enddo                 enddo
3880  c$$$               print*,'nclsx         ',nclsx  c$$$               print*,'nclsx         ',nclsx
# Line 3878  c$$$               print*,'xs(2,nclsx)   Line 3885  c$$$               print*,'xs(2,nclsx)  
3885              else                          !=== Y views              else                          !=== Y views
3886                 nclsy = nclsy + 1                 nclsy = nclsy + 1
3887                 planey(nclsy) = ip                 planey(nclsy) = ip
3888                 sgnlys(nclsy) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3889                   if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy)
3890                 clsy(nclsy)   = icl                 clsy(nclsy)   = icl
3891                 do is=1,2                 do is=1,2
3892  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
3893                    call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)  c                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
3894                      call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.,0.,0.)
3895                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2
3896                 enddo                 enddo
3897  c$$$               print*,'nclsy         ',nclsy  c$$$               print*,'nclsy         ',nclsy
# Line 3892  c$$$               print*,'ys(1,nclsy)   Line 3901  c$$$               print*,'ys(1,nclsy)  
3901  c$$$               print*,'ys(2,nclsy)   ',ys(2,nclsy)  c$$$               print*,'ys(2,nclsy)   ',ys(2,nclsy)
3902              endif              endif
3903           endif           endif
 c      print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip)  
3904    
3905  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO
3906           whichtrack(icl) = cl_used(icl)           whichtrack(icl) = cl_used(icl)
3907    *     --------------------------------------------------
3908    *     per non perdere la testa...
3909    *     whichtrack(icl) e` una variabile del common level1
3910    *     che serve solo per sapere quali cluster sono stati
3911    *     associati ad una traccia, e permettere di salvare
3912    *     solo questi nell'albero di uscita
3913    *     --------------------------------------------------
3914            
3915    
3916    c$$$         print*,' cl ',icl,' --> ',cl_used(icl)
3917    c$$$
3918    c$$$         if( cl_used(icl).ne.0 )then
3919    c$$$            if(
3920    c$$$     $           mod(VIEW(icl),2).eq.0.and.
3921    c$$$     $           cltrx(ip,whichtrack(icl)).ne.icl )
3922    c$$$     $           print*,'**WARNING** cltrx(',ip,',',whichtrack(icl)
3923    c$$$     $           ,')=',cltrx(ip,whichtrack(icl)),'.ne.',icl
3924    c$$$            if(
3925    c$$$     $           mod(VIEW(icl),2).eq.1.and.
3926    c$$$     $           cltry(ip,whichtrack(icl)).ne.icl )
3927    c$$$     $           print*,'**WARNING** cltry(',ip,',',whichtrack(icl)
3928    c$$$     $           ,')=',cltry(ip,whichtrack(icl)),'.ne.',icl
3929    c$$$         endif
3930            
3931    
3932        enddo        enddo
3933        end        end
3934    
3935    ***************************************************
3936    *                                                 *
3937    *                                                 *
3938    *                                                 *
3939    *                                                 *
3940    *                                                 *
3941    *                                                 *
3942    **************************************************
3943    
3944          subroutine fill_hough
3945    
3946    *     -------------------------------------------------------
3947    *     This routine fills the  variables related to the hough
3948    *     transform, for the debig n-tuple
3949    *     -------------------------------------------------------
3950    
3951          include 'commontracker.f'
3952          include 'level1.f'
3953          include 'common_momanhough.f'
3954          include 'common_hough.f'
3955          include 'level2.f'
3956    
3957          if(.false.
3958         $     .or.ntrpt.gt.ntrpt_max_nt
3959         $     .or.ndblt.gt.ndblt_max_nt
3960         $     .or.NCLOUDS_XZ.gt.ncloxz_max
3961         $     .or.NCLOUDS_yZ.gt.ncloyz_max
3962         $     )then
3963             ntrpt_nt=0
3964             ndblt_nt=0
3965             NCLOUDS_XZ_nt=0
3966             NCLOUDS_YZ_nt=0
3967          else
3968             ndblt_nt=ndblt
3969             ntrpt_nt=ntrpt
3970             if(ndblt.ne.0)then
3971                do id=1,ndblt
3972                   alfayz1_nt(id)=alfayz1(id) !Y0
3973                   alfayz2_nt(id)=alfayz2(id) !tg theta-yz
3974                enddo
3975             endif
3976             if(ndblt.ne.0)then
3977                do it=1,ntrpt
3978                   alfaxz1_nt(it)=alfaxz1(it) !X0
3979                   alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz
3980                   alfaxz3_nt(it)=alfaxz3(it) !1/r
3981                enddo
3982             endif
3983             nclouds_yz_nt=nclouds_yz
3984             nclouds_xz_nt=nclouds_xz
3985             if(nclouds_yz.ne.0)then
3986                nnn=0
3987                do iyz=1,nclouds_yz
3988                   ptcloud_yz_nt(iyz)=ptcloud_yz(iyz)
3989                   alfayz1_av_nt(iyz)=alfayz1_av(iyz)
3990                   alfayz2_av_nt(iyz)=alfayz2_av(iyz)
3991                   nnn=nnn+ptcloud_yz(iyz)
3992                enddo
3993                do ipt=1,nnn
3994                   db_cloud_nt(ipt)=db_cloud(ipt)
3995                 enddo
3996             endif
3997             if(nclouds_xz.ne.0)then
3998                nnn=0
3999                do ixz=1,nclouds_xz
4000                   ptcloud_xz_nt(ixz)=ptcloud_xz(ixz)
4001                   alfaxz1_av_nt(ixz)=alfaxz1_av(ixz)
4002                   alfaxz2_av_nt(ixz)=alfaxz2_av(ixz)
4003                   alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)
4004                   nnn=nnn+ptcloud_xz(ixz)              
4005                enddo
4006                do ipt=1,nnn
4007                  tr_cloud_nt(ipt)=tr_cloud(ipt)
4008                 enddo
4009             endif
4010          endif
4011          end
4012          
Legend:
Removed from v.1.8  
changed lines
  Added in v.1.31
  ViewVC Help
Powered by ViewVC 1.1.23