/[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.11 by pam-fi, Tue Nov 7 15:55:11 2006 UTC revision 1.36 by pam-fi, Tue Nov 25 14:41:38 2008 UTC
# Line 18  Line 18 
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'
 c      include 'level1.f'  
21        include 'level2.f'        include 'level2.f'
22    
23  c      include 'momanhough_init.f'  
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
# Line 42  c      include 'momanhough_init.f' Line 74  c      include 'momanhough_init.f'
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               !go to next event           goto 880               !go to next event
# Line 76  c      iflag=0 Line 108  c      iflag=0
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               !go to next event                       goto 880               !go to next event            
# Line 183  c      iflag=0 Line 215  c      iflag=0
215  c      include 'momanhough_init.f'  c      include 'momanhough_init.f'
216                
217        logical FIMAGE            !        logical FIMAGE            !
218          real trackimage(NTRACKSMAX)
219        real*8 AL_GUESS(5)        real*8 AL_GUESS(5)
220    
221  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
# Line 234  c$$$            endif Line 267  c$$$            endif
267  c$$$         enddo  c$$$         enddo
268  c$$$         if(ibest.eq.0)goto 880 !>> no good candidates  c$$$         if(ibest.eq.0)goto 880 !>> no good candidates
269    
270           rchi2best=1000000000.  *     -------------------------------------------------------
271           ndofbest=0             !(1)  *     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=0             !(1)               ndof=ndof        
281               do ii=1,nplanes    !(1)       $            +int(xgood_store(ii,i))
282                 ndof=ndof        !(1)       $            +int(ygood_store(ii,i))
283       $              +int(xgood_store(ii,i)) !(1)             enddo              
284       $              +int(ygood_store(ii,i)) !(1)             if(ndof.gt.ndofbest)then
285               enddo              !(1)               ibest=i
286               if(ndof.ge.ndofbest)then !(1)               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                 ndofbest=ndof    !(1)                 ndofbest=ndof  
294               endif              !(1)               endif            
295             endif             endif
296           enddo           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 270  c$$$         if(ibest.eq.0)goto 880 !>> Line 332  c$$$         if(ibest.eq.0)goto 880 !>>
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 286  c$$$         if(ibest.eq.0)goto 880 !>> Line 350  c$$$         if(ibest.eq.0)goto 880 !>>
350           do i=1,5           do i=1,5
351              AL_GUESS(i)=AL(i)              AL_GUESS(i)=AL(i)
352           enddo           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))            
# Line 296  c$$$         if(ibest.eq.0)goto 880 !>> Line 361  c$$$         if(ibest.eq.0)goto 880 !>>
361           jstep=0                !# minimization steps           jstep=0                !# minimization steps
362    
363           iprint=0           iprint=0
364  c         if(DEBUG)iprint=1  c         if(DEBUG.EQ.1)iprint=1
365           if(.true.)iprint=1           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(.true.)then              if(VERBOSE.EQ.1)then
370                 print *,                 print *,
371       $              '*** MINIMIZATION FAILURE *** (after refinement) '       $              '*** MINIMIZATION FAILURE *** (after refinement) '
372       $              ,iev       $              ,iev
373    
374                 print*,'guess:   ',(al_guess(i),i=1,5)  c$$$               print*,'guess:   ',(al_guess(i),i=1,5)
375                 print*,'previous: ',(al_store(i,icand),i=1,5)  c$$$               print*,'previous: ',(al_store(i,icand),i=1,5)
376                 print*,'result:   ',(al(i),i=1,5)  c$$$               print*,'result:   ',(al(i),i=1,5)
377                 print*,'xgood ',xgood  c$$$               print*,'xgood ',xgood
378                 print*,'ygood ',ygood  c$$$               print*,'ygood ',ygood
379                 print*,'----------------------------------------------'  c$$$               print*,'----------------------------------------------'
380              endif              endif
381  c            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 326  c            chi2=-chi2 Line 392  c            chi2=-chi2
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 342  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 373  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 409  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 423  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 619  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 657  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'
675        include 'level1.f'        include 'level1.f'
676        include 'calib.f'        include 'calib.f'
 c      include 'level1.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'
 c      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 xi,yi,zi
692          double precision xi_A,yi_A,zi_A
693          double precision xi_B,yi_B,zi_B
694        double precision xrt,yrt,zrt        double precision xrt,yrt,zrt
695        double precision xrt_A,yrt_A,zrt_A        double precision xrt_A,yrt_A,zrt_A
696        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  
697                
698    
699        parameter (ndivx=30)        parameter (ndivx=30)
700    
701    
702    c$$$      print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy
703                
704        resxPAM = 0        resxPAM = 0
705        resyPAM = 0        resyPAM = 0
706    
707        xPAM = 0.        xPAM = 0.D0
708        yPAM = 0.        yPAM = 0.D0
709        zPAM = 0.        zPAM = 0.D0
710        xPAM_A = 0.        xPAM_A = 0.D0
711        yPAM_A = 0.        yPAM_A = 0.D0
712        zPAM_A = 0.        zPAM_A = 0.D0
713        xPAM_B = 0.        xPAM_B = 0.D0
714        yPAM_B = 0.        yPAM_B = 0.D0
715        zPAM_B = 0.        zPAM_B = 0.D0
716    cc      print*,'## xyz_PAM: ',icx,icy,sensor,PFAx,PFAy,angx,angy
717    
718          if(sensor.lt.1.or.sensor.gt.2)then
719             print*,'xyz_PAM   ***ERROR*** wrong input '
720             print*,'sensor ',sensor
721             icx=0
722             icy=0
723          endif
724    
725  *     -----------------  *     -----------------
726  *     CLUSTER X  *     CLUSTER X
727  *     -----------------  *     -----------------      
   
728        if(icx.ne.0)then        if(icx.ne.0)then
729           viewx = VIEW(icx)  
730           nldx = nld(MAXS(icx),VIEW(icx))           viewx   = VIEW(icx)
731           nplx = npl(VIEW(icx))           nldx    = nld(MAXS(icx),VIEW(icx))
732           resxPAM = RESXAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!           nplx    = npl(VIEW(icx))
733            c         resxPAM = RESXAV
734           stripx = float(MAXS(icx))           stripx  = float(MAXS(icx))
735           if(PFAx.eq.'COG1')then  !(1)  
736              stripx = stripx      !(1)           if(
737              resxPAM = resxPAM    !(1)       $        viewx.lt.1.or.        
738           elseif(PFAx.eq.'COG2')then       $        viewx.gt.12.or.
739              stripx = stripx + cog(2,icx)                   $        nldx.lt.1.or.
740              resxPAM = resxPAM*fbad_cog(2,icx)       $        nldx.gt.3.or.
741           elseif(PFAx.eq.'ETA2')then       $        stripx.lt.1.or.
742  c            cog2 = cog(2,icx)       $        stripx.gt.3072.or.
743  c            etacorr = pfaeta2(cog2,viewx,nldx,angx)                   $        .false.)then
744  c            stripx = stripx + etacorr              print*,'xyz_PAM   ***ERROR*** wrong input '
745              stripx = stripx + pfaeta2(icx,angx)            !(3)              print*,'icx ',icx,'view ',viewx,'nld ',nldx,'strip ',stripx
746              resxPAM = risx_eta2(angx)                       !   (4)              icx = 0
747              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  
748           endif           endif
749    
750        endif  *        --------------------------
751  c      if(icy.eq.0.and.icx.ne.0)  *        magnetic-field corrections
752  c     $     print*,PFAx,icx,angx,stripx,resxPAM,'***'  *        --------------------------
753                   stripx  = stripx +  fieldcorr(viewx,bfy)      
754             angx    = effectiveangle(ax,viewx,bfy)
755            
756             call applypfa(PFAx,icx,angx,corr,res)
757             stripx  = stripx + corr
758             resxPAM = res
759    
760     10   endif
761        
762  *     -----------------  *     -----------------
763  *     CLUSTER Y  *     CLUSTER Y
764  *     -----------------  *     -----------------
765    
766        if(icy.ne.0)then        if(icy.ne.0)then
767    
768           viewy = VIEW(icy)           viewy = VIEW(icy)
769           nldy = nld(MAXS(icy),VIEW(icy))           nldy = nld(MAXS(icy),VIEW(icy))
770           nply = npl(VIEW(icy))           nply = npl(VIEW(icy))
771           resyPAM = RESYAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!  c         resyPAM = RESYAV
772             stripy = float(MAXS(icy))
773    
774             if(
775         $        viewy.lt.1.or.        
776         $        viewy.gt.12.or.
777         $        nldy.lt.1.or.
778         $        nldy.gt.3.or.
779         $        stripy.lt.1.or.
780         $        stripy.gt.3072.or.
781         $        .false.)then
782                print*,'xyz_PAM   ***ERROR*** wrong input '
783                print*,'icy ',icy,'view ',viewy,'nld ',nldy,'strip ',stripy
784                icy = 0
785                goto 20
786             endif
787    
788           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
789              print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '              if(DEBUG.EQ.1) then
790       $           ,icx,icy                 print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '
791         $              ,icx,icy
792                endif
793              goto 100              goto 100
794           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  
795    
796        endif  *        --------------------------
797    *        magnetic-field corrections
798    *        --------------------------
799             stripy  = stripy +  fieldcorr(viewy,bfx)      
800             angy    = effectiveangle(ay,viewy,bfx)
801            
802             call applypfa(PFAy,icy,angy,corr,res)
803             stripy  = stripy + corr
804             resyPAM = res
805    
806     20   endif
807    
808    cc      print*,'## stripx,stripy ',stripx,stripy
809    
         
810  c===========================================================  c===========================================================
811  C     COUPLE  C     COUPLE
812  C===========================================================  C===========================================================
# Line 840  c     (xi,yi,zi) = mechanical coordinate Line 817  c     (xi,yi,zi) = mechanical coordinate
817  c------------------------------------------------------------------------  c------------------------------------------------------------------------
818           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
819       $        .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...
820              print*,'xyz_PAM (couple):',              if(DEBUG.EQ.1) then
821       $          ' WARNING: false X strip: strip ',stripx                 print*,'xyz_PAM (couple):',
822         $              ' WARNING: false X strip: strip ',stripx
823                endif
824           endif           endif
825           xi = acoordsi(stripx,viewx)           xi = dcoordsi(stripx,viewx)
826           yi = acoordsi(stripy,viewy)           yi = dcoordsi(stripy,viewy)
827           zi = 0.           zi = 0.D0
828                    
   
829  c------------------------------------------------------------------------  c------------------------------------------------------------------------
830  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
831  c------------------------------------------------------------------------  c------------------------------------------------------------------------
# Line 881  c--------------------------------------- Line 859  c---------------------------------------
859           yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4           yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4
860           zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4           zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4
861    
862           xPAM_A = 0.           xPAM_A = 0.D0
863           yPAM_A = 0.           yPAM_A = 0.D0
864           zPAM_A = 0.           zPAM_A = 0.D0
865    
866           xPAM_B = 0.           xPAM_B = 0.D0
867           yPAM_B = 0.           yPAM_B = 0.D0
868           zPAM_B = 0.           zPAM_B = 0.D0
869    
870        elseif(        elseif(
871       $        (icx.ne.0.and.icy.eq.0).or.       $        (icx.ne.0.and.icy.eq.0).or.
# Line 907  C======================================= Line 885  C=======================================
885              nldx = nldy              nldx = nldy
886              viewx = viewy + 1              viewx = viewy + 1
887    
888              yi   = acoordsi(stripy,viewy)              yi   = dcoordsi(stripy,viewy)
889    
890              xi_A = edgeY_d - SiDimX/2              xi_A = edgeY_d - SiDimX/2
891              yi_A = yi              yi_A = yi
# Line 933  c            print*,'X-singlet ',icx,npl Line 911  c            print*,'X-singlet ',icx,npl
911  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...
912              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
913       $           .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...
914                 print*,'xyz_PAM (X-singlet):',                 if(DEBUG.EQ.1) then
915       $             ' WARNING: false X strip: strip ',stripx                    print*,'xyz_PAM (X-singlet):',
916         $                 ' WARNING: false X strip: strip ',stripx
917                   endif
918              endif              endif
919              xi   = acoordsi(stripx,viewx)              xi   = dcoordsi(stripx,viewx)
920    
921              xi_A = xi              xi_A = xi
922              yi_A = edgeX_d - SiDimY/2              yi_A = edgeX_d - SiDimY/2
# Line 956  c            print*,'X-cl ',icx,stripx,' Line 936  c            print*,'X-cl ',icx,stripx,'
936  c            print*,yi_A,' <--> ',yi_B  c            print*,yi_A,' <--> ',yi_B
937    
938           else           else
939                if(DEBUG.EQ.1) then
940              print *,'routine xyz_PAM ---> not properly used !!!'                 print *,'routine xyz_PAM ---> not properly used !!!'
941              print *,'icx = ',icx                 print *,'icx = ',icx
942              print *,'icy = ',icy                 print *,'icy = ',icy
943                endif
944              goto 100              goto 100
945                            
946           endif           endif
# Line 1008  c     (xPAM,yPAM,zPAM) = measured coordi Line 989  c     (xPAM,yPAM,zPAM) = measured coordi
989  c                        in PAMELA reference system  c                        in PAMELA reference system
990  c------------------------------------------------------------------------  c------------------------------------------------------------------------
991    
992           xPAM = 0.           xPAM = 0.D0
993           yPAM = 0.           yPAM = 0.D0
994           zPAM = 0.           zPAM = 0.D0
995    
996           xPAM_A = dcoord(xrt_A,viewx,nldx,sensor) / 1.d4           xPAM_A = dcoord(xrt_A,viewx,nldx,sensor) / 1.d4
997           yPAM_A = dcoord(yrt_A,viewy,nldy,sensor) / 1.d4           yPAM_A = dcoord(yrt_A,viewy,nldy,sensor) / 1.d4
# Line 1024  c--------------------------------------- Line 1005  c---------------------------------------
1005  c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'  c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'
1006    
1007        else        else
1008                       if(DEBUG.EQ.1) then
1009           print *,'routine xyz_PAM ---> not properly used !!!'              print *,'routine xyz_PAM ---> not properly used !!!'
1010           print *,'icx = ',icx              print *,'icx = ',icx
1011           print *,'icy = ',icy              print *,'icy = ',icy
1012                         endif
1013        endif        endif
1014                    
1015    
1016    c      print*,'## xPAM,yPAM,zPAM       ',xPAM,yPAM,zPAM
1017    c      print*,'## xPAM_A,yPAM_A,zPAM_A ',xPAM_A,yPAM_A,zPAM_A
1018    c      print*,'## xPAM_B,yPAM_B,zPAM_B ',xPAM_B,yPAM_B,zPAM_B
1019    
1020   100  continue   100  continue
1021        end        end
1022    
1023    ************************************************************************
1024    *     Call xyz_PAM subroutine with default PFA and fill the mini2 common.
1025    *     (done to be called from c/c++)
1026    ************************************************************************
1027    
1028          subroutine xyzpam(ip,icx,icy,lad,sensor,ax,ay,bfx,bfy)
1029    
1030          include 'commontracker.f'
1031          include 'level1.f'
1032          include 'common_mini_2.f'
1033          include 'common_xyzPAM.f'
1034          include 'common_mech.f'
1035          include 'calib.f'
1036          
1037    *     flag to chose PFA
1038    c$$$      character*10 PFA
1039    c$$$      common/FINALPFA/PFA
1040    
1041          integer icx,icy           !X-Y cluster ID
1042          integer sensor
1043          character*4 PFAx,PFAy     !PFA to be used
1044          real ax,ay                !X-Y geometric angle
1045          real bfx,bfy              !X-Y b-field components
1046    
1047          ipx=0
1048          ipy=0      
1049          
1050    c$$$      PFAx = 'COG4'!PFA
1051    c$$$      PFAy = 'COG4'!PFA
1052    
1053    
1054          if(icx.gt.nclstr1.or.icy.gt.nclstr1)then
1055                print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1056         $           ,' do not exists (n.clusters=',nclstr1,')'
1057                icx = -1*icx
1058                icy = -1*icy
1059                return
1060            
1061          endif
1062          
1063          call idtoc(pfaid,PFAx)
1064          call idtoc(pfaid,PFAy)
1065    
1066    cc      print*,PFAx,PFAy
1067    
1068    c$$$      call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1069    
1070    c$$$      print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy
1071          
1072          if(icx.ne.0.and.icy.ne.0)then
1073    
1074             ipx=npl(VIEW(icx))
1075             ipy=npl(VIEW(icy))
1076    c$$$         if( (nplanes-ipx+1).ne.ip.or.(nplanes-ipy+1).ne.ip )
1077    c$$$     $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1078    c$$$     $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1079            
1080             if( (nplanes-ipx+1).ne.ip )then
1081                print*,'xyzpam: ***WARNING*** cluster ',icx
1082         $           ,' does not belong to plane: ',ip
1083                icx = -1*icx
1084                return
1085             endif
1086             if( (nplanes-ipy+1).ne.ip )then
1087                print*,'xyzpam: ***WARNING*** cluster ',icy
1088         $           ,' does not belong to plane: ',ip
1089                icy = -1*icy
1090                return
1091             endif
1092    
1093             call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1094    
1095             xgood(ip) = 1.
1096             ygood(ip) = 1.
1097             resx(ip)  = resxPAM
1098             resy(ip)  = resyPAM
1099    
1100             xm(ip) = xPAM
1101             ym(ip) = yPAM
1102             zm(ip) = zPAM
1103             xm_A(ip) = 0.D0
1104             ym_A(ip) = 0.D0
1105             xm_B(ip) = 0.D0
1106             ym_B(ip) = 0.D0
1107    
1108    c         zv(ip) = zPAM
1109    
1110          elseif(icx.eq.0.and.icy.ne.0)then
1111    
1112             ipy=npl(VIEW(icy))
1113    c$$$         if((nplanes-ipy+1).ne.ip)
1114    c$$$     $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1115    c$$$     $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1116             if( (nplanes-ipy+1).ne.ip )then
1117                print*,'xyzpam: ***WARNING*** cluster ',icy
1118         $           ,' does not belong to plane: ',ip
1119                icy = -1*icy
1120                return
1121             endif
1122    
1123             call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1124            
1125             xgood(ip) = 0.
1126             ygood(ip) = 1.
1127             resx(ip)  = 1000.
1128             resy(ip)  = resyPAM
1129    
1130    cPP --- old ---
1131    c$$$         xm(ip) = -100.
1132    c$$$         ym(ip) = -100.
1133    c$$$         zm(ip) = (zPAM_A+zPAM_B)/2.
1134    c$$$         xm_A(ip) = xPAM_A
1135    c$$$         ym_A(ip) = yPAM_A
1136    c$$$         xm_B(ip) = xPAM_B
1137    c$$$         ym_B(ip) = yPAM_B
1138    cPP --- new ---
1139             xm(ip) = -100.
1140             ym(ip) = -100.
1141             zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1142             xm_A(ip) = xPAM_A
1143             ym_A(ip) = yPAM_A
1144             zm_A(ip) = zPAM_A
1145             xm_B(ip) = xPAM_B
1146             ym_B(ip) = yPAM_B
1147             zm_B(ip) = zPAM_B
1148    cPP -----------
1149    
1150    c         zv(ip) = (zPAM_A+zPAM_B)/2.
1151            
1152          elseif(icx.ne.0.and.icy.eq.0)then
1153    
1154             ipx=npl(VIEW(icx))
1155    c$$$         if((nplanes-ipx+1).ne.ip)
1156    c$$$     $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1157    c$$$     $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1158    
1159             if( (nplanes-ipx+1).ne.ip )then
1160                print*,'xyzpam: ***WARNING*** cluster ',icx
1161         $           ,' does not belong to plane: ',ip
1162                icx = -1*icx
1163                return
1164             endif
1165    
1166             call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1167          
1168             xgood(ip) = 1.
1169             ygood(ip) = 0.
1170             resx(ip)  = resxPAM
1171             resy(ip)  = 1000.
1172    
1173    cPP --- old ---
1174    c$$$         xm(ip) = -100.
1175    c$$$         ym(ip) = -100.
1176    c$$$         zm(ip) = (zPAM_A+zPAM_B)/2.
1177    c$$$         xm_A(ip) = xPAM_A
1178    c$$$         ym_A(ip) = yPAM_A
1179    c$$$         xm_B(ip) = xPAM_B
1180    c$$$         ym_B(ip) = yPAM_B
1181    cPP --- new ---
1182             xm(ip) = -100.
1183             ym(ip) = -100.
1184             zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1185             xm_A(ip) = xPAM_A
1186             ym_A(ip) = yPAM_A
1187             zm_A(ip) = zPAM_A
1188             xm_B(ip) = xPAM_B
1189             ym_B(ip) = yPAM_B
1190             zm_B(ip) = zPAM_B
1191    cPP -----------
1192    
1193    c         zv(ip) = (zPAM_A+zPAM_B)/2.
1194    
1195          else
1196    
1197             il = 2
1198             if(lad.ne.0)il=lad
1199             is = 1
1200             if(sensor.ne.0)is=sensor
1201    c         print*,nplanes-ip+1,il,is
1202    
1203             xgood(ip) = 0.
1204             ygood(ip) = 0.
1205             resx(ip)  = 1000.
1206             resy(ip)  = 1000.
1207    
1208             xm(ip) = -100.
1209             ym(ip) = -100.          
1210             zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1211             xm_A(ip) = 0.
1212             ym_A(ip) = 0.
1213             xm_B(ip) = 0.
1214             ym_B(ip) = 0.
1215    
1216    c         zv(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1217    
1218          endif
1219    
1220          if(DEBUG.EQ.1)then
1221    c         print*,'----------------------------- track coord'
1222    22222    format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
1223             write(*,22222)ip,zm(ip),xm(ip),ym(ip)
1224         $        ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip)
1225         $        ,xgood(ip),ygood(ip),resx(ip),resy(ip)
1226    c$$$         print*,'-----------------------------'
1227          endif
1228          end
1229    
1230  ********************************************************************************  ********************************************************************************
1231  ********************************************************************************  ********************************************************************************
# Line 1055  c         print*,'A-(',xPAM_A,yPAM_A,') Line 1247  c         print*,'A-(',xPAM_A,yPAM_A,')
1247  *      *    
1248  ********************************************************************************  ********************************************************************************
1249    
1250        real function distance_to(XPP,YPP)        real function distance_to(rXPP,rYPP)
1251    
1252        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
1253    
# Line 1064  c         print*,'A-(',xPAM_A,yPAM_A,') Line 1256  c         print*,'A-(',xPAM_A,yPAM_A,')
1256  *     ( i.e. distance/resolution )  *     ( i.e. distance/resolution )
1257  *     -----------------------------------  *     -----------------------------------
1258    
1259          real rXPP,rYPP
1260          double precision XPP,YPP
1261        double precision distance,RE        double precision distance,RE
1262        double precision BETA,ALFA,xmi,ymi        double precision BETA,ALFA,xmi,ymi
1263    
1264          XPP=DBLE(rXPP)
1265          YPP=DBLE(rYPP)
1266    
1267  *     ----------------------  *     ----------------------
1268        if (        if (
1269       +     xPAM.eq.0.and.       +     xPAM.eq.0.and.
# Line 1109  c         print*,'A-(',xPAM_A,yPAM_A,') Line 1306  c         print*,'A-(',xPAM_A,yPAM_A,')
1306           endif                   endif        
1307    
1308           distance=           distance=
1309       $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2       $       ((xmi-XPP)**2+(ymi-YPP)**2)!QUIQUI
1310    cc     $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2
1311           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1312    
1313  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 1134  c$$$         print*,' resolution ',re Line 1332  c$$$         print*,' resolution ',re
1332  *     ----------------------  *     ----------------------
1333                    
1334           distance=           distance=
1335       $        ((xPAM-XPP)/resxPAM)**2       $       ((xPAM-XPP))**2 !QUIQUI
1336       $        +       $       +
1337       $        ((yPAM-YPP)/resyPAM)**2       $       ((yPAM-YPP))**2
1338    c$$$     $        ((xPAM-XPP)/resxPAM)**2
1339    c$$$     $        +
1340    c$$$     $        ((yPAM-YPP)/resyPAM)**2
1341           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1342    
1343  c$$$         print*,xPAM,yPAM,zPAM  c$$$         print*,xPAM,yPAM,zPAM
# Line 1145  c$$$         print*,' resolution ',resxP Line 1346  c$$$         print*,' resolution ',resxP
1346                    
1347        else        else
1348                    
1349           print*  c         print*
1350       $        ,' function distance_to ---> wrong usage!!!'  c     $        ,' function distance_to ---> wrong usage!!!'
1351           print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM  c         print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM
1352           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 '
1353       $        ,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
1354        endif          endif  
1355    
1356        distance_to = sngl(distance)        distance_to = sngl(distance)
# Line 1189  c$$$         print*,' resolution ',resxP Line 1390  c$$$         print*,' resolution ',resxP
1390        data c1/1.,0.,0.,1./        data c1/1.,0.,0.,1./
1391        data c2/1.,-1.,-1.,1./        data c2/1.,-1.,-1.,1./
1392        data c3/1.,1.,0.,0./        data c3/1.,1.,0.,0./
1393        real*8 yvvv,xvvv        double precision yvvv,xvvv
1394        double precision xi,yi,zi        double precision xi,yi,zi
1395        double precision xrt,yrt,zrt        double precision xrt,yrt,zrt
1396        real AA,BB        real AA,BB
1397        real yvv(4),xvv(4)        double precision yvv(4),xvv(4)
1398    
1399  *     tollerance to consider the track inside the sensitive area  *     tollerance to consider the track inside the sensitive area
1400        real ptoll        real ptoll
1401        data ptoll/150./          !um        data ptoll/150./          !um
1402    
1403        external nviewx,nviewy,acoordsi,dcoord        external nviewx,nviewy,dcoordsi,dcoord
1404    
1405        nplpt = nplPAM            !plane        nplpt = nplPAM            !plane
1406        viewx = nviewx(nplpt)        viewx = nviewx(nplpt)
# Line 1214  c$$$         print*,' resolution ',resxP Line 1415  c$$$         print*,' resolution ',resxP
1415  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1416  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
1417  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1418                 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)  c               if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
1419       $              .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...  c     $              .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
1420  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...
1421                    print*,'whichsensor: ',  c                  print*,'whichsensor: ',
1422       $                ' WARNING: false X strip: strip ',stripx  c     $                ' WARNING: false X strip: strip ',stripx
1423                 endif  c               endif
1424                 xi = acoordsi(stripx,viewx)                 xi = dcoordsi(stripx,viewx)
1425                 yi = acoordsi(stripy,viewy)                 yi = dcoordsi(stripy,viewy)
1426                 zi = 0.                 zi = 0.D0
1427  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1428  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
1429  c------------------------------------------------------------------------  c------------------------------------------------------------------------
# Line 1373  c      include 'common_analysis.f' Line 1574  c      include 'common_analysis.f'
1574        is_cp=0        is_cp=0
1575        if(id.lt.0)is_cp=1        if(id.lt.0)is_cp=1
1576        if(id.gt.0)is_cp=2        if(id.gt.0)is_cp=2
1577        if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'  c      if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'
1578    
1579        return        return
1580        end        end
# Line 1444  c      include 'common_analysis.f' Line 1645  c      include 'common_analysis.f'
1645  *************************************************************************  *************************************************************************
1646  *************************************************************************  *************************************************************************
1647  *************************************************************************  *************************************************************************
 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  
1648                
1649    
1650  ***************************************************  ***************************************************
# Line 1740  c      include 'level1.f' Line 1673  c      include 'level1.f'
1673        integer iflag        integer iflag
1674    
1675        integer badseed,badclx,badcly        integer badseed,badclx,badcly
1676        
1677          if(DEBUG.EQ.1)print*,'cl_to_couples:'
1678    
1679  *     init variables  *     init variables
1680        ncp_tot=0        ncp_tot=0
# Line 1769  c      include 'level1.f' Line 1704  c      include 'level1.f'
1704        enddo        enddo
1705  *     mask views with too many clusters  *     mask views with too many clusters
1706        do iv=1,nviews        do iv=1,nviews
1707           if( ncl_view(iv).gt. nclustermax)then           if( ncl_view(iv).gt. nclusterlimit)then
1708              mask_view(iv) = 1  c            mask_view(iv) = 1
1709              if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '              mask_view(iv) = mask_view(iv) + 2**0
1710       $           ,nclustermax,' on view ', iv,' --> masked!'              if(DEBUG.EQ.1)
1711         $        print*,' * WARNING * cl_to_couple: n.clusters > '
1712         $        ,nclusterlimit,' on view ', iv,' --> masked!'
1713           endif           endif
1714        enddo        enddo
1715    
# Line 1789  c      include 'level1.f' Line 1726  c      include 'level1.f'
1726  *     ----------------------------------------------------  *     ----------------------------------------------------
1727  *     cut on charge (X VIEW)  *     cut on charge (X VIEW)
1728  *     ----------------------------------------------------  *     ----------------------------------------------------
1729           if(dedx(icx).lt.dedx_x_min)then           if(sgnl(icx).lt.dedx_x_min)then
1730              cl_single(icx)=0              cl_single(icx)=0
1731              goto 10              goto 10
1732           endif           endif
# Line 1839  c     endif Line 1776  c     endif
1776  *     ----------------------------------------------------  *     ----------------------------------------------------
1777  *     cut on charge (Y VIEW)  *     cut on charge (Y VIEW)
1778  *     ----------------------------------------------------  *     ----------------------------------------------------
1779              if(dedx(icy).lt.dedx_y_min)then              if(sgnl(icy).lt.dedx_y_min)then
1780                 cl_single(icy)=0                 cl_single(icy)=0
1781                 goto 20                 goto 20
1782              endif              endif
# Line 1885  c     endif Line 1822  c     endif
1822  *     charge correlation  *     charge correlation
1823  *     (modified to be applied only below saturation... obviously)  *     (modified to be applied only below saturation... obviously)
1824    
1825                 if(  .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx)                 if(  .not.(sgnl(icy).gt.chsaty.and.sgnl(icx).gt.chsatx)
1826       $              .and.       $              .and.
1827       $              .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx)       $              .not.(sgnl(icy).lt.chmipy.and.sgnl(icx).lt.chmipx)
1828       $              .and.       $              .and.
1829       $              (badclx.eq.1.and.badcly.eq.1)       $              (badclx.eq.1.and.badcly.eq.1)
1830       $              .and.       $              .and.
1831       $              .true.)then       $              .true.)then
1832    
1833                    ddd=(dedx(icy)                    ddd=(sgnl(icy)
1834       $                 -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))       $                 -kch(nplx,nldx)*sgnl(icx)-cch(nplx,nldx))
1835                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
1836    
1837  c                  cut = chcut * sch(nplx,nldx)  c                  cut = chcut * sch(nplx,nldx)
1838    
1839                    sss=(kch(nplx,nldx)*dedx(icy)+dedx(icx)                    sss=(kch(nplx,nldx)*sgnl(icy)+sgnl(icx)
1840       $                 -kch(nplx,nldx)*cch(nplx,nldx))       $                 -kch(nplx,nldx)*cch(nplx,nldx))
1841                    sss=sss/sqrt(kch(nplx,nldx)**2+1)                    sss=sss/sqrt(kch(nplx,nldx)**2+1)
1842                    cut = chcut * (16 + sss/50.)                    cut = chcut * (16 + sss/50.)
# Line 1908  c                  cut = chcut * sch(npl Line 1845  c                  cut = chcut * sch(npl
1845                       goto 20    !charge not consistent                       goto 20    !charge not consistent
1846                    endif                    endif
1847                 endif                 endif
1848                  
1849  *     ------------------> COUPLE <------------------                 if(ncp_plane(nplx).gt.ncouplemax)then
1850  *     check to do not overflow vector dimentions                    if(verbose.eq.1)print*,
 c$$$               if(ncp_plane(nplx).gt.ncouplemax)then  
 c$$$                  if(DEBUG)print*,  
 c$$$     $                    ' ** warning ** number of identified'//  
 c$$$     $                    ' couples on plane ',nplx,  
 c$$$     $                    ' exceeds vector dimention'//  
 c$$$     $                    ' ( ',ncouplemax,' )'  
 c$$$c     good2=.false.  
 c$$$c     goto 880   !fill ntp and go to next event  
 c$$$                  iflag=1  
 c$$$                  return  
 c$$$               endif  
                 
                if(ncp_plane(nplx).eq.ncouplemax)then  
                   if(verbose)print*,  
1851       $                 '** warning ** number of identified '//       $                 '** warning ** number of identified '//
1852       $                 'couples on plane ',nplx,       $                 'couples on plane ',nplx,
1853       $                 'exceeds vector dimention '       $                 'exceeds vector dimention '
1854       $                 ,'( ',ncouplemax,' )'       $                 ,'( ',ncouplemax,' ) --> masked!'
1855  c     good2=.false.  c                  mask_view(nviewx(nplx)) = 2
1856  c     goto 880   !fill ntp and go to next event  c                  mask_view(nviewy(nply)) = 2
1857                    mask_view(nviewx(nplx)) = 2                    mask_view(nviewx(nplx))= mask_view(nviewx(nplx))+ 2**1
1858                    mask_view(nviewy(nply)) = 2                    mask_view(nviewy(nply))= mask_view(nviewy(nply))+ 2**1
1859  c                  iflag=1                    goto 10
 c                  return  
1860                 endif                 endif
1861                                
1862    *     ------------------> COUPLE <------------------
1863                 ncp_plane(nplx) = ncp_plane(nplx) + 1                 ncp_plane(nplx) = ncp_plane(nplx) + 1
1864                 clx(nplx,ncp_plane(nplx))=icx                 clx(nplx,ncp_plane(nplx))=icx
1865                 cly(nply,ncp_plane(nplx))=icy                 cly(nply,ncp_plane(nplx))=icy
1866                 cl_single(icx)=0                 cl_single(icx)=0
1867                 cl_single(icy)=0                 cl_single(icy)=0
             endif                                
1868  *     ----------------------------------------------  *     ----------------------------------------------
1869    
1870                endif                              
1871    
1872   20         continue   20         continue
1873           enddo                  !end loop on clusters(Y)           enddo                  !end loop on clusters(Y)
1874                    
# Line 1961  c                  return Line 1885  c                  return
1885        enddo        enddo
1886                
1887                
1888        if(DEBUG)then        if(DEBUG.EQ.1)then
1889           print*,'clusters  ',nclstr1           print*,'clusters  ',nclstr1
1890           print*,'good    ',(cl_good(i),i=1,nclstr1)           print*,'good    ',(cl_good(i),i=1,nclstr1)
1891           print*,'singles ',(cl_single(i),i=1,nclstr1)           print*,'singlets',(cl_single(i),i=1,nclstr1)
1892           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
1893        endif        endif
1894                
1895        do ip=1,6        do ip=1,6
1896           ncp_tot = ncp_tot + ncp_plane(ip)           ncp_tot = ncp_tot + ncp_plane(ip)
1897        enddo        enddo
 c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  
         
 c$$$      if(ncp_tot.gt.ncp_max)then  
 c$$$         if(verbose)print*,  
 c$$$     $           '** warning ** number of identified '//  
 c$$$     $           'couples exceeds upper limit for Hough tr. '  
 c$$$     $           ,'( ',ncp_max,' )'              
 c$$$         iflag=1  
 c$$$         return  
 c$$$      endif  
1898                
1899        return        return
1900        end        end
# Line 1993  c$$$      endif Line 1907  c$$$      endif
1907  *                                                 *  *                                                 *
1908  *                                                 *  *                                                 *
1909  **************************************************  **************************************************
 c$$$      subroutine cl_to_couples_nocharge(iflag)  
 c$$$  
 c$$$      include 'commontracker.f'  
 c$$$      include 'level1.f'  
 c$$$      include 'common_momanhough.f'  
 c$$$c      include 'momanhough_init.f'  
 c$$$      include 'calib.f'  
 c$$$c      include 'level1.f'  
 c$$$  
 c$$$  
 c$$$*     output flag  
 c$$$*     --------------  
 c$$$*     0 = good event  
 c$$$*     1 = bad event  
 c$$$*     --------------  
 c$$$      integer iflag  
 c$$$  
 c$$$      integer badseed,badcl  
 c$$$  
 c$$$*     init variables  
 c$$$      ncp_tot=0  
 c$$$      do ip=1,nplanes  
 c$$$         do ico=1,ncouplemax  
 c$$$            clx(ip,ico)=0  
 c$$$            cly(ip,ico)=0  
 c$$$         enddo  
 c$$$         ncp_plane(ip)=0  
 c$$$         do icl=1,nclstrmax_level2  
 c$$$            cls(ip,icl)=1  
 c$$$         enddo  
 c$$$         ncls(ip)=0  
 c$$$      enddo  
 c$$$      do icl=1,nclstrmax_level2  
 c$$$         cl_single(icl)=1  
 c$$$         cl_good(icl)=0  
 c$$$      enddo  
 c$$$        
 c$$$*     start association  
 c$$$      ncouples=0  
 c$$$      do icx=1,nclstr1          !loop on cluster (X)  
 c$$$         if(mod(VIEW(icx),2).eq.1)goto 10  
 c$$$          
 c$$$*     ----------------------------------------------------  
 c$$$*     cut on charge (X VIEW)  
 c$$$         if(dedx(icx).lt.dedx_x_min)then  
 c$$$            cl_single(icx)=0  
 c$$$            goto 10  
 c$$$         endif  
 c$$$*     cut BAD (X VIEW)              
 c$$$         badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))  
 c$$$         ifirst=INDSTART(icx)  
 c$$$         if(icx.ne.nclstr1) then  
 c$$$            ilast=INDSTART(icx+1)-1  
 c$$$         else  
 c$$$            ilast=TOTCLLENGTH  
 c$$$         endif  
 c$$$         badcl=badseed  
 c$$$         do igood=-ngoodstr,ngoodstr  
 c$$$            ibad=1  
 c$$$            if((INDMAX(icx)+igood).gt.ifirst.and.  
 c$$$     $           (INDMAX(icx)+igood).lt.ilast.and.  
 c$$$     $           .true.)then  
 c$$$               ibad=BAD(VIEW(icx),  
 c$$$     $              nvk(MAXS(icx)+igood),  
 c$$$     $              nst(MAXS(icx)+igood))  
 c$$$            endif  
 c$$$            badcl=badcl*ibad  
 c$$$         enddo  
 c$$$         if(badcl.eq.0)then     !<<<<<<<<<<<<<< BAD cut  
 c$$$            cl_single(icx)=0    !<<<<<<<<<<<<<< BAD cut  
 c$$$            goto 10             !<<<<<<<<<<<<<< BAD cut  
 c$$$         endif                  !<<<<<<<<<<<<<< BAD cut  
 c$$$*     ----------------------------------------------------  
 c$$$          
 c$$$         cl_good(icx)=1  
 c$$$         nplx=npl(VIEW(icx))  
 c$$$         nldx=nld(MAXS(icx),VIEW(icx))  
 c$$$          
 c$$$         do icy=1,nclstr1       !loop on cluster (Y)  
 c$$$            if(mod(VIEW(icy),2).eq.0)goto 20  
 c$$$              
 c$$$*     ----------------------------------------------------  
 c$$$*     cut on charge (Y VIEW)  
 c$$$            if(dedx(icy).lt.dedx_y_min)then  
 c$$$               cl_single(icy)=0  
 c$$$               goto 20  
 c$$$            endif  
 c$$$*     cut BAD (Y VIEW)              
 c$$$            badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))  
 c$$$            ifirst=INDSTART(icy)  
 c$$$            if(icy.ne.nclstr1) then  
 c$$$               ilast=INDSTART(icy+1)-1  
 c$$$            else  
 c$$$               ilast=TOTCLLENGTH  
 c$$$            endif  
 c$$$            badcl=badseed  
 c$$$            do igood=-ngoodstr,ngoodstr  
 c$$$               ibad=1  
 c$$$               if((INDMAX(icy)+igood).gt.ifirst.and.  
 c$$$     $              (INDMAX(icy)+igood).lt.ilast.and.  
 c$$$     $              .true.)  
 c$$$     $              ibad=BAD(VIEW(icy),  
 c$$$     $              nvk(MAXS(icy)+igood),  
 c$$$     $              nst(MAXS(icy)+igood))  
 c$$$               badcl=badcl*ibad  
 c$$$            enddo  
 c$$$            if(badcl.eq.0)then  !<<<<<<<<<<<<<< BAD cut  
 c$$$               cl_single(icy)=0 !<<<<<<<<<<<<<< BAD cut  
 c$$$               goto 20          !<<<<<<<<<<<<<< BAD cut  
 c$$$            endif               !<<<<<<<<<<<<<< BAD cut  
 c$$$*     ----------------------------------------------------  
 c$$$              
 c$$$              
 c$$$            cl_good(icy)=1                    
 c$$$            nply=npl(VIEW(icy))  
 c$$$            nldy=nld(MAXS(icy),VIEW(icy))  
 c$$$              
 c$$$*     ----------------------------------------------  
 c$$$*     CONDITION TO FORM A COUPLE  
 c$$$*     ----------------------------------------------  
 c$$$*     geometrical consistency (same plane and ladder)  
 c$$$            if(nply.eq.nplx.and.nldy.eq.nldx)then  
 c$$$*     charge correlation  
 c$$$*     ===========================================================  
 c$$$*     this version of the subroutine is used for the calibration  
 c$$$*     thus charge-correlation selection is obviously removed  
 c$$$*     ===========================================================  
 c$$$c$$$               ddd=(dedx(icy)  
 c$$$c$$$     $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))  
 c$$$c$$$               ddd=ddd/sqrt(kch(nplx,nldx)**2+1)  
 c$$$c$$$               cut=chcut*sch(nplx,nldx)  
 c$$$c$$$               if(abs(ddd).gt.cut)goto 20 !charge not consistent  
 c$$$*     ===========================================================  
 c$$$                
 c$$$                
 c$$$*     ------------------> COUPLE <------------------  
 c$$$*     check to do not overflow vector dimentions  
 c$$$c$$$               if(ncp_plane(nplx).gt.ncouplemax)then  
 c$$$c$$$                  if(DEBUG)print*,  
 c$$$c$$$     $                    ' ** warning ** number of identified'//  
 c$$$c$$$     $                    ' couples on plane ',nplx,  
 c$$$c$$$     $                    ' exceeds vector dimention'//  
 c$$$c$$$     $                    ' ( ',ncouplemax,' )'  
 c$$$c$$$c     good2=.false.  
 c$$$c$$$c     goto 880   !fill ntp and go to next event  
 c$$$c$$$                  iflag=1  
 c$$$c$$$                  return  
 c$$$c$$$               endif  
 c$$$                
 c$$$               if(ncp_plane(nplx).eq.ncouplemax)then  
 c$$$                  if(verbose)print*,  
 c$$$     $                 '** warning ** number of identified '//  
 c$$$     $                 'couples on plane ',nplx,  
 c$$$     $                 'exceeds vector dimention '  
 c$$$     $                 ,'( ',ncouplemax,' )'  
 c$$$c     good2=.false.  
 c$$$c     goto 880   !fill ntp and go to next event                      
 c$$$                  iflag=1  
 c$$$                  return  
 c$$$               endif  
 c$$$                
 c$$$               ncp_plane(nplx) = ncp_plane(nplx) + 1  
 c$$$               clx(nplx,ncp_plane(nplx))=icx  
 c$$$               cly(nply,ncp_plane(nplx))=icy  
 c$$$               cl_single(icx)=0  
 c$$$               cl_single(icy)=0  
 c$$$            endif                                
 c$$$*     ----------------------------------------------  
 c$$$  
 c$$$ 20         continue  
 c$$$         enddo                  !end loop on clusters(Y)  
 c$$$          
 c$$$ 10      continue  
 c$$$      enddo                     !end loop on clusters(X)  
 c$$$        
 c$$$        
 c$$$      do icl=1,nclstr1  
 c$$$         if(cl_single(icl).eq.1)then  
 c$$$            ip=npl(VIEW(icl))  
 c$$$            ncls(ip)=ncls(ip)+1  
 c$$$            cls(ip,ncls(ip))=icl  
 c$$$         endif  
 c$$$      enddo  
 c$$$        
 c$$$        
 c$$$      if(DEBUG)then  
 c$$$         print*,'clusters  ',nclstr1  
 c$$$         print*,'good    ',(cl_good(i),i=1,nclstr1)  
 c$$$         print*,'singles ',(cl_single(i),i=1,nclstr1)  
 c$$$         print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)  
 c$$$      endif  
 c$$$        
 c$$$      do ip=1,6  
 c$$$         ncp_tot=ncp_tot+ncp_plane(ip)  
 c$$$      enddo  
 c$$$c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  
 c$$$        
 c$$$      if(ncp_tot.gt.ncp_max)then  
 c$$$         if(verbose)print*,  
 c$$$     $           '** warning ** number of identified '//  
 c$$$     $           'couples exceeds upper limit for Hough tr. '  
 c$$$     $           ,'( ',ncp_max,' )'              
 c$$$c            good2=.false.  
 c$$$c     goto 880       !fill ntp and go to next event  
 c$$$         iflag=1  
 c$$$         return  
 c$$$      endif  
 c$$$        
 c$$$      return  
 c$$$      end  
 c$$$  
         
 ***************************************************  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 **************************************************  
1910    
1911        subroutine cp_to_doubtrip(iflag)        subroutine cp_to_doubtrip(iflag)
 c*****************************************************  
 c     02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
1912    
1913        include 'commontracker.f'        include 'commontracker.f'
1914        include 'level1.f'        include 'level1.f'
1915        include 'common_momanhough.f'        include 'common_momanhough.f'
 c      include 'momanhough_init.f'  
1916        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
1917        include 'common_mini_2.f'        include 'common_mini_2.f'
1918        include 'calib.f'        include 'calib.f'
 c      include 'level1.f'  
1919    
1920    
1921  *     output flag  *     output flag
# Line 2259  c      double precision xm3,ym3,zm3 Line 1948  c      double precision xm3,ym3,zm3
1948        real xc,zc,radius        real xc,zc,radius
1949  *     -----------------------------  *     -----------------------------
1950    
1951          if(DEBUG.EQ.1)print*,'cp_to_doubtrip:'
1952    
1953    *     --------------------------------------------
1954    *     put a limit to the maximum number of couples
1955    *     per plane, in order to apply hough transform
1956    *     (couples recovered during track refinement)
1957    *     --------------------------------------------
1958          do ip=1,nplanes
1959             if(ncp_plane(ip).gt.ncouplelimit)then
1960    c            mask_view(nviewx(ip)) = 8
1961    c            mask_view(nviewy(ip)) = 8
1962                mask_view(nviewx(ip)) = mask_view(nviewx(ip)) + 2**7
1963                mask_view(nviewy(ip)) = mask_view(nviewy(ip)) + 2**7
1964             endif
1965          enddo
1966    
1967    
1968        ndblt=0                   !number of doublets        ndblt=0                   !number of doublets
1969        ntrpt=0                   !number of triplets        ntrpt=0                   !number of triplets
1970                
1971        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1
1972           do is1=1,2             !loop on sensors - COPPIA 1           if(  mask_view(nviewx(ip1)).ne.0 .or.
1973                     $        mask_view(nviewy(ip1)).ne.0 )goto 10 !skip plane
1974             do is1=1,2             !loop on sensors - COPPIA 1            
1975              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1
1976                 icx1=clx(ip1,icp1)                 icx1=clx(ip1,icp1)
1977                 icy1=cly(ip1,icp1)                 icy1=cly(ip1,icp1)
1978  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)
1979                 call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)  c               call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)
1980                   call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.)
1981                 xm1=xPAM                 xm1=xPAM
1982                 ym1=yPAM                 ym1=yPAM
1983                 zm1=zPAM                                   zm1=zPAM                  
1984  c     print*,'***',is1,xm1,ym1,zm1  c     print*,'***',is1,xm1,ym1,zm1
1985    
1986                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
1987                    do is2=1,2    !loop on sensors -ndblt COPPIA 2                    if(  mask_view(nviewx(ip2)).ne.0 .or.
1988                             $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane
1989                      do is2=1,2    !loop on sensors -ndblt COPPIA 2                    
1990                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
1991                          icx2=clx(ip2,icp2)                          icx2=clx(ip2,icp2)
1992                          icy2=cly(ip2,icp2)                          icy2=cly(ip2,icp2)
1993  c                        call xyz_PAM  c                        call xyz_PAM
1994  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)
1995    c                        call xyz_PAM
1996    c     $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)
1997                          call xyz_PAM                          call xyz_PAM
1998       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.)
1999                          xm2=xPAM                          xm2=xPAM
2000                          ym2=yPAM                          ym2=yPAM
2001                          zm2=zPAM                          zm2=zPAM
# Line 2295  c     $                       (icx2,icy2 Line 2005  c     $                       (icx2,icy2
2005  *     (2 couples needed)  *     (2 couples needed)
2006  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2007                          if(ndblt.eq.ndblt_max)then                          if(ndblt.eq.ndblt_max)then
2008                             if(verbose)print*,                             if(verbose.eq.1)print*,
2009       $                          '** warning ** number of identified '//       $                          '** warning ** number of identified '//
2010       $                          'doublets exceeds vector dimention '       $                          'doublets exceeds vector dimention '
2011       $                          ,'( ',ndblt_max,' )'       $                          ,'( ',ndblt_max,' )'
2012  c                           good2=.false.  c                           good2=.false.
2013  c                           goto 880 !fill ntp and go to next event  c                           goto 880 !fill ntp and go to next event
2014                             do iv=1,12                             do iv=1,12
2015                                mask_view(iv) = 3  c                              mask_view(iv) = 3
2016                                  mask_view(iv) = mask_view(iv)+ 2**2
2017                             enddo                             enddo
2018                             iflag=1                             iflag=1
2019                             return                             return
# Line 2337  c$$$ Line 2048  c$$$
2048  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2049    
2050    
2051                          if(ip2.eq.nplanes)goto 30 !no possible combination with 3 couples                          if(ip2.eq.nplanes)goto 31 !no possible combination with 3 couples
2052    
2053                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3
2054                               if(  mask_view(nviewx(ip3)).ne.0 .or.
2055         $                          mask_view(nviewy(ip3)).ne.0 )goto 30 !skip plane
2056                             do is3=1,2 !loop on sensors - COPPIA 3                             do is3=1,2 !loop on sensors - COPPIA 3
2057                                                                
2058                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3
# Line 2346  c$$$ Line 2060  c$$$
2060                                   icy3=cly(ip3,icp3)                                   icy3=cly(ip3,icp3)
2061  c                                 call xyz_PAM  c                                 call xyz_PAM
2062  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)
2063    c                                 call xyz_PAM
2064    c     $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)
2065                                   call xyz_PAM                                   call xyz_PAM
2066       $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)       $                                (icx3,icy3,is3,PFAdef,PFAdef
2067         $                                ,0.,0.,0.,0.)
2068                                   xm3=xPAM                                   xm3=xPAM
2069                                   ym3=yPAM                                   ym3=yPAM
2070                                   zm3=zPAM                                   zm3=zPAM
2071  *     find the circle passing through the three points  *     find the circle passing through the three points
2072    c$$$                                 call tricircle(3,xp,zp,angp,resp,chi
2073    c$$$     $                                ,xc,zc,radius,iflag)
2074                                     iflag = DEBUG
2075                                   call tricircle(3,xp,zp,angp,resp,chi                                   call tricircle(3,xp,zp,angp,resp,chi
2076       $                                ,xc,zc,radius,iflag)       $                                ,xc,zc,radius,iflag)
2077  c     print*,xc,zc,radius  c     print*,xc,zc,radius
# Line 2368  c     $                                 Line 2088  c     $                                
2088  *     (3 couples needed)  *     (3 couples needed)
2089  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2090                                   if(ntrpt.eq.ntrpt_max)then                                   if(ntrpt.eq.ntrpt_max)then
2091                                      if(verbose)print*,                                      if(verbose.eq.1)print*,
2092       $                     '** warning ** number of identified '//       $                     '** warning ** number of identified '//
2093       $                     'triplets exceeds vector dimention '       $                     'triplets exceeds vector dimention '
2094       $                    ,'( ',ntrpt_max,' )'       $                    ,'( ',ntrpt_max,' )'
2095  c                                    good2=.false.  c                                    good2=.false.
2096  c                                    goto 880 !fill ntp and go to next event  c                                    goto 880 !fill ntp and go to next event
2097                                      do iv=1,nviews                                      do iv=1,nviews
2098                                         mask_view(iv) = 4  c                                       mask_view(iv) = 4
2099                                           mask_view(iv)=mask_view(iv)+ 2**3
2100                                      enddo                                      enddo
2101                                      iflag=1                                      iflag=1
2102                                      return                                      return
# Line 2415  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp Line 2136  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp
2136                                endif                                endif
2137                             enddo !end loop on COPPIA 3                             enddo !end loop on COPPIA 3
2138                          enddo   !end loop on sensors - COPPIA 3                          enddo   !end loop on sensors - COPPIA 3
2139     30                     continue
2140                       enddo      !end loop on planes  - COPPIA 3                       enddo      !end loop on planes  - COPPIA 3
2141   30                  continue   31                  continue
2142                                            
2143   1                enddo         !end loop on COPPIA 2   1                enddo         !end loop on COPPIA 2
2144                 enddo            !end loop on sensors - COPPIA 2                 enddo            !end loop on sensors - COPPIA 2
2145     20            continue
2146              enddo               !end loop on planes  - COPPIA 2              enddo               !end loop on planes  - COPPIA 2
2147                            
2148           enddo                  !end loop on COPPIA1           enddo                  !end loop on COPPIA1
2149        enddo                     !end loop on sensors - COPPIA 1        enddo                     !end loop on sensors - COPPIA 1
2150     10   continue
2151        enddo                     !end loop on planes  - COPPIA 1        enddo                     !end loop on planes  - COPPIA 1
2152                
2153        if(DEBUG)then        if(DEBUG.EQ.1)then
2154           print*,'--- doublets ',ndblt           print*,'--- doublets ',ndblt
2155           print*,'--- triplets ',ntrpt           print*,'--- triplets ',ntrpt
2156        endif        endif
# Line 2473  c      include 'momanhough_init.f' Line 2197  c      include 'momanhough_init.f'
2197        integer cp_useds1(ncouplemaxtot) ! sensor 1        integer cp_useds1(ncouplemaxtot) ! sensor 1
2198        integer cp_useds2(ncouplemaxtot) ! sensor 2        integer cp_useds2(ncouplemaxtot) ! sensor 2
2199    
2200          if(DEBUG.EQ.1)print*,'doub_to_YZcloud:'
2201    
2202  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2203  *     classification of DOUBLETS  *     classification of DOUBLETS
# Line 2599  c         if(ncpused.lt.ncpyz_min)goto 2 Line 2324  c         if(ncpused.lt.ncpyz_min)goto 2
2324  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2325    
2326           if(nclouds_yz.ge.ncloyz_max)then           if(nclouds_yz.ge.ncloyz_max)then
2327              if(verbose)print*,              if(verbose.eq.1)print*,
2328       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2329       $           'YZ clouds exceeds vector dimention '       $           'YZ clouds exceeds vector dimention '
2330       $           ,'( ',ncloyz_max,' )'       $           ,'( ',ncloyz_max,' )'
2331  c               good2=.false.  c               good2=.false.
2332  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2333              do iv=1,nviews              do iv=1,nviews
2334                 mask_view(iv) = 5  c               mask_view(iv) = 5
2335                   mask_view(iv) = mask_view(iv) + 2**4
2336              enddo              enddo
2337              iflag=1              iflag=1
2338              return              return
# Line 2626  c     ptcloud_yz_nt(nclouds_yz)=npt Line 2352  c     ptcloud_yz_nt(nclouds_yz)=npt
2352  c     print*,'>> ',ipt,db_all(ipt)  c     print*,'>> ',ipt,db_all(ipt)
2353           enddo             enddo  
2354           npt_tot=npt_tot+npt           npt_tot=npt_tot+npt
2355           if(DEBUG)then           if(DEBUG.EQ.1)then
2356              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
2357              print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'              print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'
2358              print*,'- alfayz1 ',alfayz1_av(nclouds_yz)              print*,'- alfayz1  ',alfayz1_av(nclouds_yz)
2359              print*,'- alfayz2 ',alfayz2_av(nclouds_yz)              print*,'- alfayz2  ',alfayz2_av(nclouds_yz)
2360              print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)              print*,'cp_useds1  ',(cp_useds1(icp),icp=1,ncp_tot)
2361              print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)              print*,'cp_useds2  ',(cp_useds2(icp),icp=1,ncp_tot)
2362              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)              print*,'cpcloud_yz '
2363         $           ,(cpcloud_yz(nclouds_yz,icp),icp=1,ncp_tot)
2364                print*,'hit_plane  ',(hit_plane(ip),ip=1,nplanes)
2365  c$$$            print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = '  c$$$            print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = '
2366  c$$$     $           ,ptcloud_yz(nclouds_yz)  c$$$     $           ,ptcloud_yz(nclouds_yz)
2367  c$$$            print*,'nt-uple: db_cloud(...) = '  c$$$            print*,'nt-uple: db_cloud(...) = '
# Line 2651  c$$$     $           ,(db_cloud(iii),iii Line 2379  c$$$     $           ,(db_cloud(iii),iii
2379          goto 90                          goto 90                
2380        endif                            endif                    
2381                
2382        if(DEBUG)then        if(DEBUG.EQ.1)then
2383           print*,'---------------------- '           print*,'---------------------- '
2384           print*,'Y-Z total clouds ',nclouds_yz           print*,'Y-Z total clouds ',nclouds_yz
2385           print*,' '           print*,' '
# Line 2700  c      include 'momanhough_init.f' Line 2428  c      include 'momanhough_init.f'
2428        integer cp_useds1(ncouplemaxtot) ! sensor 1        integer cp_useds1(ncouplemaxtot) ! sensor 1
2429        integer cp_useds2(ncouplemaxtot) ! sensor 2        integer cp_useds2(ncouplemaxtot) ! sensor 2
2430    
2431          if(DEBUG.EQ.1)print*,'trip_to_XZcloud:'
2432    
2433  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2434  *     classification of TRIPLETS  *     classification of TRIPLETS
2435  *     according to distance in parameter space  *     according to distance in parameter space
# Line 2757  c         tr_temp(1)=itr1 Line 2487  c         tr_temp(1)=itr1
2487       $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2                     $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2              
2488                 distance = sqrt(distance)                 distance = sqrt(distance)
2489                                
2490                 if(distance.lt.cutdistxz)then  *     ------------------------------------------------------------------------
2491    *     FORCE INCLUSION OF TRIPLETS COMPOSED BY SAME COUPLES, IGNORING THE IMAGE
2492    *     ------------------------------------------------------------------------
2493    *     (added in august 2007)
2494                   istrimage=0
2495                   if(
2496         $              abs(cpxz1(itrref)).eq.abs(cpxz1(itr2)).and.
2497         $              abs(cpxz2(itrref)).eq.abs(cpxz2(itr2)).and.
2498         $              abs(cpxz3(itrref)).eq.abs(cpxz3(itr2)).and.
2499         $              .true.)istrimage=1
2500    
2501                   if(distance.lt.cutdistxz.or.istrimage.eq.1)then
2502  c     print*,idb1,idb2,distance,' cloud ',nclouds_yz  c     print*,idb1,idb2,distance,' cloud ',nclouds_yz
2503                    if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1                    if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1
2504                    if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1                    if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1
# Line 2816  c     print*,'check cp_used' Line 2557  c     print*,'check cp_used'
2557           enddo           enddo
2558  c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet  c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet
2559           if(npt.lt.nptxz_min)goto 22288     !next triplet           if(npt.lt.nptxz_min)goto 22288     !next triplet
2560           if(nplused.lt.nplxz_min)goto 22288 !next doublet           if(nplused.lt.nplxz_min)goto 22288 !next triplet
2561                    
2562  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2563  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2564           if(nclouds_xz.ge.ncloxz_max)then           if(nclouds_xz.ge.ncloxz_max)then
2565              if(verbose)print*,              if(verbose.eq.1)print*,
2566       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2567       $           'XZ clouds exceeds vector dimention '       $           'XZ clouds exceeds vector dimention '
2568       $           ,'( ',ncloxz_max,' )'       $           ,'( ',ncloxz_max,' )'
2569  c     good2=.false.  c     good2=.false.
2570  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2571              do iv=1,nviews              do iv=1,nviews
2572                 mask_view(iv) = 6  c               mask_view(iv) = 6
2573                   mask_view(iv) =  mask_view(iv) + 2**5
2574              enddo              enddo
2575              iflag=1              iflag=1
2576              return              return
# Line 2847  c     goto 880         !fill ntp and go Line 2589  c     goto 880         !fill ntp and go
2589           enddo           enddo
2590           npt_tot=npt_tot+npt           npt_tot=npt_tot+npt
2591                    
2592           if(DEBUG)then           if(DEBUG.EQ.1)then
2593              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
2594              print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'                            print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'              
2595              print*,'- alfaxz1 ',alfaxz1_av(nclouds_xz)              print*,'- alfaxz1  ',alfaxz1_av(nclouds_xz)
2596              print*,'- alfaxz2 ',alfaxz2_av(nclouds_xz)              print*,'- alfaxz2  ',alfaxz2_av(nclouds_xz)
2597              print*,'- alfaxz3 ',alfaxz3_av(nclouds_xz)              print*,'- alfaxz3  ',alfaxz3_av(nclouds_xz)
2598              print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)              print*,'cp_useds1  ',(cp_useds1(icp),icp=1,ncp_tot)
2599              print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)              print*,'cp_useds2  ',(cp_useds2(icp),icp=1,ncp_tot)
2600                print*,'cpcloud_xz '
2601         $           ,(cpcloud_xz(nclouds_xz,icp),icp=1,ncp_tot)
2602              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)
2603  c$$$            print*,'nt-uple: ptcloud_xz(',nclouds_xz,') = '  c$$$            print*,'nt-uple: ptcloud_xz(',nclouds_xz,') = '
2604  c$$$     $           ,ptcloud_xz(nclouds_xz)  c$$$     $           ,ptcloud_xz(nclouds_xz)
# Line 2872  c$$$     $           ,(tr_cloud(iii),iii Line 2616  c$$$     $           ,(tr_cloud(iii),iii
2616           goto 91                           goto 91                
2617         endif                             endif                    
2618                
2619        if(DEBUG)then        if(DEBUG.EQ.1)then
2620           print*,'---------------------- '           print*,'---------------------- '
2621           print*,'X-Z total clouds ',nclouds_xz           print*,'X-Z total clouds ',nclouds_xz
2622           print*,' '           print*,' '
# Line 2893  c$$$     $           ,(tr_cloud(iii),iii Line 2637  c$$$     $           ,(tr_cloud(iii),iii
2637  **************************************************  **************************************************
2638    
2639        subroutine clouds_to_ctrack(iflag)        subroutine clouds_to_ctrack(iflag)
 c*****************************************************  
 c     02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
2640    
2641        include 'commontracker.f'        include 'commontracker.f'
2642        include 'level1.f'        include 'level1.f'
# Line 2903  c*************************************** Line 2644  c***************************************
2644        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2645        include 'common_mini_2.f'        include 'common_mini_2.f'
2646        include 'common_mech.f'        include 'common_mech.f'
2647  c      include 'momanhough_init.f'  
2648    
2649    
2650  *     output flag  *     output flag
# Line 2927  c      include 'momanhough_init.f' Line 2668  c      include 'momanhough_init.f'
2668  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2669  *     variables for track fitting  *     variables for track fitting
2670        double precision AL_INI(5)        double precision AL_INI(5)
 c      double precision tath  
2671  *     -----------------------------------------------------------  *     -----------------------------------------------------------
 c      real fitz(nplanes)        !z coordinates of the planes in cm  
2672    
2673          if(DEBUG.EQ.1)print*,'clouds_to_ctrack:'
2674    
2675    
2676        ntracks=0                 !counter of track candidates        ntracks=0                 !counter of track candidates
# Line 2952  c      real fitz(nplanes)        !z coor Line 2692  c      real fitz(nplanes)        !z coor
2692                 enddo                 enddo
2693              enddo              enddo
2694              ncp_ok=0              ncp_ok=0
2695              do icp=1,ncp_tot    !loop on couples              do icp=1,ncp_tot    !loop over couples
2696  *     get info on  *     get info on
2697                 cpintersec(icp)=min(                 cpintersec(icp)=min(
2698       $              cpcloud_yz(iyz,icp),       $              cpcloud_yz(iyz,icp),
# Line 2961  c      real fitz(nplanes)        !z coor Line 2701  c      real fitz(nplanes)        !z coor
2701       $    (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.
2702       $    (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.
2703       $              .false.)cpintersec(icp)=0       $              .false.)cpintersec(icp)=0
2704    *     cpintersec is >0 if yz and xz clouds contain the same image of couple icp
2705                 if(cpintersec(icp).ne.0)then                 if(cpintersec(icp).ne.0)then
2706                    ncp_ok=ncp_ok+1                      ncp_ok=ncp_ok+1  
2707                                        
# Line 2993  c      real fitz(nplanes)        !z coor Line 2734  c      real fitz(nplanes)        !z coor
2734                 nplused=nplused+ hit_plane(ip)                 nplused=nplused+ hit_plane(ip)
2735              enddo              enddo
2736                            
             if(nplused.lt.nplxz_min)goto 888 !next doublet  
             if(ncp_ok.lt.ncpok)goto 888 !next cloud  
2737                            
2738              if(DEBUG)then              if(DEBUG.EQ.1)then
2739                 print*,'Combination ',iyz,ixz                 print*,'Combination ',iyz,ixz
2740       $              ,' db ',ptcloud_yz(iyz)       $              ,' db ',ptcloud_yz(iyz)
2741       $              ,' tr ',ptcloud_xz(ixz)       $              ,' tr ',ptcloud_xz(ixz)
2742       $              ,'  -----> # matching couples ',ncp_ok       $              ,'  -----> # matching couples ',ncp_ok
2743              endif              endif
2744    
2745    c            if(nplused.lt.nplxz_min)goto 888 !next combination
2746                if(nplused.lt.nplyz_min)goto 888 !next combination
2747                if(ncp_ok.lt.ncpok)goto 888 !next combination
2748    
2749  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
2750  c$$$  print*,'Configurazione cluster XZ'  c$$$  print*,'Configurazione cluster XZ'
2751  c$$$  print*,'1 -- ',(clx(1,i),i=1,ncp_plane(1))  c$$$  print*,'1 -- ',(clx(1,i),i=1,ncp_plane(1))
# Line 3047  c$$$            AL_INI(5) = (1.e2*alfaxz Line 2791  c$$$            AL_INI(5) = (1.e2*alfaxz
2791  c$$$              c$$$            
2792  c$$$            if(AL_INI(5).gt.defmax)goto 888 !next cloud  c$$$            if(AL_INI(5).gt.defmax)goto 888 !next cloud
2793                                                    
2794              if(DEBUG)then              if(DEBUG.EQ.1)then
2795                   print*,'track candidate', ntracks+1
2796                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))
2797                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))
2798                 print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4))                 print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4))
# Line 3080  c$$$            if(AL_INI(5).gt.defmax)g Line 2825  c$$$            if(AL_INI(5).gt.defmax)g
2825                                hit_plane(6)=icp6                                hit_plane(6)=icp6
2826                                if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6                                if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6
2827                                                                
2828                                  *                             ---------------------------------------
2829    *                             check if this group of couples has been
2830    *                             already fitted    
2831    *                             ---------------------------------------
2832                                  do ica=1,ntracks
2833                                     isthesame=1
2834                                     do ip=1,NPLANES
2835                                        if(hit_plane(ip).ne.0)then
2836                                           if(  CP_STORE(nplanes-ip+1,ica)
2837         $                                      .ne.
2838         $                                      cp_match(ip,hit_plane(ip)) )
2839         $                                      isthesame=0
2840                                        else
2841                                           if(  CP_STORE(nplanes-ip+1,ica)
2842         $                                      .ne.
2843         $                                      0 )
2844         $                                      isthesame=0
2845                                        endif
2846                                     enddo
2847                                     if(isthesame.eq.1)then
2848                                        if(DEBUG.eq.1)
2849         $                                   print*,'(already fitted)'
2850                                        goto 666 !jump to next combination
2851                                     endif
2852                                  enddo
2853    
2854                                call track_init !init TRACK common                                call track_init !init TRACK common
2855    
2856                                do ip=1,nplanes !loop on planes                                do ip=1,nplanes !loop on planes (bottom to top)
2857                                   if(hit_plane(ip).ne.0)then                                   if(hit_plane(ip).ne.0)then
2858                                      id=cp_match(ip,hit_plane(ip))                                      id=cp_match(ip,hit_plane(ip))
2859                                      is=is_cp(id)                                      is=is_cp(id)
# Line 3097  c$$$            if(AL_INI(5).gt.defmax)g Line 2867  c$$$            if(AL_INI(5).gt.defmax)g
2867  *                                   *************************  *                                   *************************
2868  c                                    call xyz_PAM(icx,icy,is,  c                                    call xyz_PAM(icx,icy,is,
2869  c     $                                   'COG2','COG2',0.,0.)  c     $                                   'COG2','COG2',0.,0.)
2870    c                                    call xyz_PAM(icx,icy,is, !(1)
2871    c     $                                   PFAdef,PFAdef,0.,0.) !(1)
2872                                      call xyz_PAM(icx,icy,is, !(1)                                      call xyz_PAM(icx,icy,is, !(1)
2873       $                                   PFAdef,PFAdef,0.,0.) !(1)       $                                   PFAdef,PFAdef,0.,0.,0.,0.)
2874  *                                   *************************  *                                   *************************
2875  *                                   -----------------------------  *                                   -----------------------------
2876                                      xgood(nplanes-ip+1)=1.                                      xgood(nplanes-ip+1)=1.
# Line 3125  c$$$                              enddo Line 2897  c$$$                              enddo
2897                                ifail=0 !error flag in chi^2 computation                                ifail=0 !error flag in chi^2 computation
2898                                jstep=0 !number of  minimization steps                                jstep=0 !number of  minimization steps
2899                                iprint=0                                iprint=0
2900  c                              if(DEBUG)iprint=1  c                              if(DEBUG.EQ.1)iprint=1
2901                                if(DEBUG)iprint=1                                if(DEBUG.EQ.1)iprint=2
2902                                call mini2(jstep,ifail,iprint)                                call mini2(jstep,ifail,iprint)
2903                                if(ifail.ne.0) then                                if(ifail.ne.0) then
2904                                   if(DEBUG)then                                   if(DEBUG.EQ.1)then
2905                                      print *,                                      print *,
2906       $                              '*** MINIMIZATION FAILURE *** '       $                              '*** MINIMIZATION FAILURE *** '
2907       $                              //'(clouds_to_ctrack)'       $                              //'(clouds_to_ctrack)'
# Line 3154  c                                 chi2=- Line 2926  c                                 chi2=-
2926  *     --------------------------  *     --------------------------
2927                                if(ntracks.eq.NTRACKSMAX)then                                if(ntracks.eq.NTRACKSMAX)then
2928                                                                    
2929                                   if(verbose)print*,                                   if(verbose.eq.1)print*,
2930       $                 '** warning ** number of candidate tracks '//       $                 '** warning ** number of candidate tracks '//
2931       $                 ' exceeds vector dimension '       $                 ' exceeds vector dimension '
2932       $                ,'( ',NTRACKSMAX,' )'       $                ,'( ',NTRACKSMAX,' )'
2933  c                                 good2=.false.  c                                 good2=.false.
2934  c                                 goto 880 !fill ntp and go to next event                      c                                 goto 880 !fill ntp and go to next event                    
2935                                   do iv=1,nviews                                   do iv=1,nviews
2936                                      mask_view(iv) = 7  c                                    mask_view(iv) = 7
2937                                        mask_view(iv) = mask_view(iv) + 2**6
2938                                   enddo                                   enddo
2939                                   iflag=1                                   iflag=1
2940                                   return                                   return
# Line 3169  c                                 goto 8 Line 2942  c                                 goto 8
2942                                                                
2943                                ntracks = ntracks + 1                                ntracks = ntracks + 1
2944                                                                
2945  c$$$                              ndof=0                                                                do ip=1,nplanes !top to bottom
2946                                do ip=1,nplanes  
 c$$$                                 ndof=ndof  
 c$$$     $                                +int(xgood(ip))  
 c$$$     $                                +int(ygood(ip))  
2947                                   XV_STORE(ip,ntracks)=sngl(xv(ip))                                   XV_STORE(ip,ntracks)=sngl(xv(ip))
2948                                   YV_STORE(ip,ntracks)=sngl(yv(ip))                                   YV_STORE(ip,ntracks)=sngl(yv(ip))
2949                                   ZV_STORE(ip,ntracks)=sngl(zv(ip))                                                                       ZV_STORE(ip,ntracks)=sngl(zv(ip))                                    
# Line 3189  c$$$     $                               Line 2959  c$$$     $                              
2959                                   AYV_STORE(ip,ntracks)=sngl(ayv(ip))                                   AYV_STORE(ip,ntracks)=sngl(ayv(ip))
2960                                   XGOOD_STORE(ip,ntracks)=sngl(xgood(ip))                                   XGOOD_STORE(ip,ntracks)=sngl(xgood(ip))
2961                                   YGOOD_STORE(ip,ntracks)=sngl(ygood(ip))                                   YGOOD_STORE(ip,ntracks)=sngl(ygood(ip))
2962    *                                NB! hit_plane is defined from bottom to top
2963                                   if(hit_plane(ip).ne.0)then                                   if(hit_plane(ip).ne.0)then
2964                                      CP_STORE(nplanes-ip+1,ntracks)=                                      CP_STORE(nplanes-ip+1,ntracks)=
2965       $                                   cp_match(ip,hit_plane(ip))       $                                   cp_match(ip,hit_plane(ip))
2966                                        SENSOR_STORE(nplanes-ip+1,ntracks)
2967         $                              = is_cp(cp_match(ip,hit_plane(ip)))
2968                                        LADDER_STORE(nplanes-ip+1,ntracks)
2969         $                                   = LADDER(
2970         $                                   clx(ip,icp_cp(
2971         $                                   cp_match(ip,hit_plane(ip)
2972         $                                   ))));
2973                                   else                                   else
2974                                      CP_STORE(nplanes-ip+1,ntracks)=0                                      CP_STORE(nplanes-ip+1,ntracks)=0
2975                                        SENSOR_STORE(nplanes-ip+1,ntracks)=0
2976                                        LADDER_STORE(nplanes-ip+1,ntracks)=0
2977                                   endif                                   endif
2978                                   CLS_STORE(nplanes-ip+1,ntracks)=0                                   BX_STORE(ip,ntracks)=0!I dont need it now
2979                                     BY_STORE(ip,ntracks)=0!I dont need it now
2980                                     CLS_STORE(ip,ntracks)=0
2981                                   do i=1,5                                   do i=1,5
2982                                      AL_STORE(i,ntracks)=sngl(AL(i))                                      AL_STORE(i,ntracks)=sngl(AL(i))
2983                                   enddo                                   enddo
2984                                enddo                                enddo
2985                                                                
 c$$$  *                             Number of Degree Of Freedom  
 c$$$  ndof=ndof-5                            
 c$$$  *                             reduced chi^2  
 c$$$  rchi2=chi2/dble(ndof)  
2986                                RCHI2_STORE(ntracks)=chi2                                RCHI2_STORE(ntracks)=chi2
2987                                                                
2988  *     --------------------------------  *     --------------------------------
# Line 3228  c$$$  rchi2=chi2/dble(ndof) Line 3006  c$$$  rchi2=chi2/dble(ndof)
3006           return           return
3007        endif        endif
3008                
3009        if(DEBUG)then  c$$$      if(DEBUG.EQ.1)then
3010           print*,'****** TRACK CANDIDATES ***********'  c$$$         print*,'****** TRACK CANDIDATES ***********'
3011           print*,'#         R. chi2        RIG'  c$$$         print*,'#         R. chi2        RIG'
3012           do i=1,ntracks  c$$$         do i=1,ntracks
3013              print*,i,' --- ',rchi2_store(i),' --- '  c$$$            print*,i,' --- ',rchi2_store(i),' --- '
3014       $           ,1./abs(AL_STORE(5,i))  c$$$     $           ,1./abs(AL_STORE(5,i))
3015           enddo  c$$$         enddo
3016           print*,'***********************************'  c$$$         print*,'***********************************'
3017    c$$$      endif
3018          if(DEBUG.EQ.1)then
3019            print*,'****** TRACK CANDIDATES *****************'
3020            print*,'#         R. chi2        RIG         ndof'
3021            do i=1,ntracks
3022              ndof=0                !(1)
3023              do ii=1,nplanes       !(1)
3024                ndof=ndof           !(1)
3025         $           +int(xgood_store(ii,i)) !(1)
3026         $           +int(ygood_store(ii,i)) !(1)
3027              enddo                 !(1)
3028              print*,i,' --- ',rchi2_store(i),' --- '
3029         $         ,1./abs(AL_STORE(5,i)),' --- ',ndof
3030            enddo
3031            print*,'*****************************************'
3032        endif        endif
3033                
3034                
# Line 3254  c$$$  rchi2=chi2/dble(ndof) Line 3047  c$$$  rchi2=chi2/dble(ndof)
3047    
3048        subroutine refine_track(ibest)        subroutine refine_track(ibest)
3049    
 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******************************************************  
3050    
3051        include 'commontracker.f'        include 'commontracker.f'
3052        include 'level1.f'        include 'level1.f'
# Line 3266  c*************************************** Line 3054  c***************************************
3054        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3055        include 'common_mini_2.f'        include 'common_mini_2.f'
3056        include 'common_mech.f'        include 'common_mech.f'
 c      include 'momanhough_init.f'  
 c      include 'level1.f'  
3057        include 'calib.f'        include 'calib.f'
3058    
   
3059  *     flag to chose PFA  *     flag to chose PFA
3060        character*10 PFA        character*10 PFA
3061        common/FINALPFA/PFA        common/FINALPFA/PFA
3062    
3063          real k(6)
3064          DATA k/1.099730,0.418900,0.220939,0.220907,0.418771,1.100674/
3065    
3066          real xp,yp,zp
3067          real xyzp(3),bxyz(3)
3068          equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3))
3069    
3070          if(DEBUG.EQ.1)print*,'refine_track:'
3071  *     =================================================  *     =================================================
3072  *     new estimate of positions using ETA algorithm  *     new estimate of positions using ETA algorithm
3073  *                          and  *                          and
# Line 3283  c      include 'level1.f' Line 3076  c      include 'level1.f'
3076        call track_init        call track_init
3077        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3078    
3079             if(DEBUG.EQ.1)print*,' ........... plane ',ip,' ........... '
3080    
3081             xP=XV_STORE(nplanes-ip+1,ibest)
3082             yP=YV_STORE(nplanes-ip+1,ibest)
3083             zP=ZV_STORE(nplanes-ip+1,ibest)
3084             call gufld(xyzp,bxyz)
3085             BX_STORE(nplanes-ip+1,ibest)=bxyz(1)
3086             BY_STORE(nplanes-ip+1,ibest)=bxyz(2)
3087    c$$$  bxyz(1)=0
3088    c$$$         bxyz(2)=0
3089    c$$$         bxyz(3)=0
3090  *     |||||||||||||||||||||||||||||||||||||||||||||||||  *     |||||||||||||||||||||||||||||||||||||||||||||||||
3091  *     -------------------------------------------------  *     -------------------------------------------------
3092  *     If the plane has been already included, it just  *     If the plane has been already included, it just
# Line 3303  c      include 'level1.f' Line 3107  c      include 'level1.f'
3107       $           ,ip_cp(id),ip       $           ,ip_cp(id),ip
3108              icx=clx(ip,icp)              icx=clx(ip,icp)
3109              icy=cly(ip,icp)              icy=cly(ip,icp)
3110    c            call xyz_PAM(icx,icy,is,
3111    c     $           PFA,PFA,
3112    c     $           AXV_STORE(nplanes-ip+1,ibest),
3113    c     $           AYV_STORE(nplanes-ip+1,ibest))
3114              call xyz_PAM(icx,icy,is,              call xyz_PAM(icx,icy,is,
 c     $           'ETA2','ETA2',  
3115       $           PFA,PFA,       $           PFA,PFA,
3116       $           AXV_STORE(nplanes-ip+1,ibest),       $           AXV_STORE(nplanes-ip+1,ibest),
3117       $           AYV_STORE(nplanes-ip+1,ibest))       $           AYV_STORE(nplanes-ip+1,ibest),
3118  c$$$  call xyz_PAM(icx,icy,is,       $           bxyz(1),
3119  c$$$  $              'COG2','COG2',       $           bxyz(2)
3120  c$$$  $              0.,       $           )
3121  c$$$  $              0.)  
3122              xm(nplanes-ip+1) = xPAM              xm(nplanes-ip+1) = xPAM
3123              ym(nplanes-ip+1) = yPAM              ym(nplanes-ip+1) = yPAM
3124              zm(nplanes-ip+1) = zPAM              zm(nplanes-ip+1) = zPAM
# Line 3320  c$$$  $              0.) Line 3127  c$$$  $              0.)
3127              resx(nplanes-ip+1) = resxPAM              resx(nplanes-ip+1) = resxPAM
3128              resy(nplanes-ip+1) = resyPAM              resy(nplanes-ip+1) = resyPAM
3129    
3130  c            dedxtrk(nplanes-ip+1) = (dedx(icx)+dedx(icy))/2. !(1)              dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx))
3131              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)  
3132                            
3133  *     |||||||||||||||||||||||||||||||||||||||||||||||||  *     |||||||||||||||||||||||||||||||||||||||||||||||||
3134  *     -------------------------------------------------  *     -------------------------------------------------
# Line 3337  c            dedxtrk(nplanes-ip+1) = (de Line 3143  c            dedxtrk(nplanes-ip+1) = (de
3143                                
3144  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3145  *     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)  
3146              call whichsensor(ip,xP,yP,nldt,ist)              call whichsensor(ip,xP,yP,nldt,ist)
3147  *     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
3148              if(nldt.eq.0.or.ist.eq.0)goto 133              if(nldt.eq.0.or.ist.eq.0)goto 133
3149    
3150                SENSOR_STORE(nplanes-ip+1,IBEST)=ist
3151                LADDER_STORE(nplanes-ip+1,IBEST)=nldt
3152  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3153    
3154              if(DEBUG)then              if(DEBUG.EQ.1)then
3155                 print*,                 print*,
3156       $              '------ Plane ',ip,' intersected on LADDER ',nldt       $              '------ Plane ',ip,' intersected on LADDER ',nldt
3157       $              ,' SENSOR ',ist       $              ,' SENSOR ',ist
# Line 3356  c            dedxtrk(nplanes-ip+1) = (de Line 3162  c            dedxtrk(nplanes-ip+1) = (de
3162  *     ===========================================  *     ===========================================
3163  *     STEP 1 >>>>>>>  try to include a new couple  *     STEP 1 >>>>>>>  try to include a new couple
3164  *     ===========================================  *     ===========================================
3165  c            if(DEBUG)print*,'>>>> try to include a new couple'  c            if(DEBUG.EQ.1)print*,'>>>> try to include a new couple'
3166              distmin=1000000.              distmin=1000000.
3167              xmm = 0.              xmm = 0.
3168              ymm = 0.              ymm = 0.
# Line 3373  c            if(DEBUG)print*,'>>>> try t Line 3179  c            if(DEBUG)print*,'>>>> try t
3179                 if(LADDER(icx).ne.nldt.or. !If the ladder number does not match                 if(LADDER(icx).ne.nldt.or. !If the ladder number does not match
3180  c     $              cl_used(icx).eq.1.or. !or the X cluster is already used  c     $              cl_used(icx).eq.1.or. !or the X cluster is already used
3181  c     $              cl_used(icy).eq.1.or. !or the Y cluster is already used  c     $              cl_used(icy).eq.1.or. !or the Y cluster is already used
3182       $              cl_used(icx).ne.0.or. !or the X cluster is already used !(3)       $              cl_used(icx).ne.0.or. !or the X cluster is already used
3183       $              cl_used(icy).ne.0.or. !or the Y cluster is already used !(3)       $              cl_used(icy).ne.0.or. !or the Y cluster is already used
3184       $              .false.)goto 1188 !then jump to next couple.       $              .false.)goto 1188 !then jump to next couple.
3185  *            *          
3186                 call xyz_PAM(icx,icy,ist,                 call xyz_PAM(icx,icy,ist,
3187       $              PFA,PFA,       $              PFA,PFA,
 c     $              'ETA2','ETA2',  
3188       $              AXV_STORE(nplanes-ip+1,ibest),       $              AXV_STORE(nplanes-ip+1,ibest),
3189       $              AYV_STORE(nplanes-ip+1,ibest))       $              AYV_STORE(nplanes-ip+1,ibest),
3190         $              bxyz(1),
3191         $              bxyz(2)
3192         $              )
3193                                
3194                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3195                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3196                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
3197                 if(DEBUG)print*,'( couple ',id                 if(DEBUG.EQ.1)print*,'( couple ',id
3198       $              ,' ) normalized distance ',distance       $              ,' ) distance ',distance
3199                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3200                    xmm = xPAM                    xmm = xPAM
3201                    ymm = yPAM                    ymm = yPAM
# Line 3396  c     $              'ETA2','ETA2', Line 3204  c     $              'ETA2','ETA2',
3204                    rymm = resyPAM                    rymm = resyPAM
3205                    distmin = distance                    distmin = distance
3206                    idm = id                                      idm = id                  
3207  c                 dedxmm = (dedx(icx)+dedx(icy))/2. !(1)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3208                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3209                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)  c     QUIQUI --> non devo moltiplicare per clinc?!?!?!
3210                      clincnewc=10*sqrt(rymm**2+rxmm**2 !QUIQUI
3211         $                 +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2))) !QUIQUI
3212                 endif                 endif
3213   1188          continue   1188          continue
3214              enddo               !end loop on couples on plane icp              enddo               !end loop on couples on plane icp
3215              if(distmin.le.clinc)then                    c            if(distmin.le.clinc)then     !QUIQUI              
3216                if(distmin.le.clincnewc)then     !QUIQUI              
3217  *              -----------------------------------  *              -----------------------------------
3218                 xm(nplanes-ip+1) = xmm         !<<<                 xm(nplanes-ip+1) = xmm !<<<
3219                 ym(nplanes-ip+1) = ymm         !<<<                 ym(nplanes-ip+1) = ymm !<<<
3220                 zm(nplanes-ip+1) = zmm         !<<<                 zm(nplanes-ip+1) = zmm !<<<
3221                 xgood(nplanes-ip+1) = 1        !<<<                 xgood(nplanes-ip+1) = 1 !<<<
3222                 ygood(nplanes-ip+1) = 1        !<<<                 ygood(nplanes-ip+1) = 1 !<<<
3223                 resx(nplanes-ip+1)=rxmm        !<<<                 resx(nplanes-ip+1)=rxmm !<<<
3224                 resy(nplanes-ip+1)=rymm        !<<<                 resy(nplanes-ip+1)=rymm !<<<
3225  c              dedxtrk(nplanes-ip+1) = dedxmm !<<<  !(1)                 dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3226                 dedxtrk_x(nplanes-ip+1) = dedxmmx    !(1)                 dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
                dedxtrk_y(nplanes-ip+1) = dedxmmy    !(1)  
3227  *              -----------------------------------  *              -----------------------------------
3228                 CP_STORE(nplanes-ip+1,ibest)=idm                       CP_STORE(nplanes-ip+1,ibest)=idm      
3229                 if(DEBUG)print*,'%%%% included couple ',idm                 if(DEBUG.EQ.1)print*,'%%%% included couple ',idm
3230       $              ,' (norm.dist.= ',distmin,', cut ',clinc,' )'       $              ,' (dist.= ',distmin,', cut ',clinc,' )'
3231                 goto 133         !next plane                 goto 133         !next plane
3232              endif              endif
3233  *     ================================================  *     ================================================
3234  *     STEP 2 >>>>>>>  try to include a single cluster  *     STEP 2 >>>>>>>  try to include a single cluster
3235  *                     either from a couple or single  *                     either from a couple or single
3236  *     ================================================  *     ================================================
3237  c            if(DEBUG)print*,'>>>> try to include a new cluster'  c            if(DEBUG.EQ.1)print*,'>>>> try to include a new cluster'
3238              distmin=1000000.              distmin=1000000.
3239              xmm_A = 0.          !---------------------------              xmm_A = 0.          !---------------------------
3240              ymm_A = 0.          ! init variables that              ymm_A = 0.          ! init variables that
# Line 3450  c            if(DEBUG)print*,'>>>> try t Line 3260  c            if(DEBUG)print*,'>>>> try t
3260  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
3261                 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)
3262  *                                              !jump to the Y cluster  *                                              !jump to the Y cluster
3263    c               call xyz_PAM(icx,0,ist,
3264    c     $              PFA,PFA,
3265    c     $              AXV_STORE(nplanes-ip+1,ibest),0.)              
3266                 call xyz_PAM(icx,0,ist,                 call xyz_PAM(icx,0,ist,
 c     $              'ETA2','ETA2',  
3267       $              PFA,PFA,       $              PFA,PFA,
3268       $              AXV_STORE(nplanes-ip+1,ibest),0.)                     $              AXV_STORE(nplanes-ip+1,ibest),0.,
3269         $              bxyz(1),
3270         $              bxyz(2)
3271         $              )              
3272                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3273                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3274                 if(DEBUG)print*,'( cl-X ',icx                 if(DEBUG.EQ.1)print*,'( cl-X ',icx
3275       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3276                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3277                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3278                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 3469  c     $              'ETA2','ETA2', Line 3284  c     $              'ETA2','ETA2',
3284                    rymm = resyPAM                    rymm = resyPAM
3285                    distmin = distance                    distmin = distance
3286                    iclm = icx                    iclm = icx
3287  c                  dedxmm = dedx(icx) !(1)  c                  dedxmm = sgnl(icx) !(1)
3288                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3289                    dedxmmy = 0.        !(1)                    dedxmmy = 0.        !(1)
3290                 endif                                   endif                  
3291  11881          continue  11881          continue
# Line 3478  c                  dedxmm = dedx(icx) !( Line 3293  c                  dedxmm = dedx(icx) !(
3293  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
3294                 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)
3295  *                                              !jump to the next couple  *                                              !jump to the next couple
3296    c               call xyz_PAM(0,icy,ist,
3297    c     $              PFA,PFA,
3298    c     $              0.,AYV_STORE(nplanes-ip+1,ibest))
3299                 call xyz_PAM(0,icy,ist,                 call xyz_PAM(0,icy,ist,
 c     $              'ETA2','ETA2',  
3300       $              PFA,PFA,       $              PFA,PFA,
3301       $              0.,AYV_STORE(nplanes-ip+1,ibest))       $              0.,AYV_STORE(nplanes-ip+1,ibest),
3302         $              bxyz(1),
3303         $              bxyz(2)
3304         $              )
3305                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3306                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3307                 if(DEBUG)print*,'( cl-Y ',icy                 if(DEBUG.EQ.1)print*,'( cl-Y ',icy
3308       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3309                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3310                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3311                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 3497  c     $              'ETA2','ETA2', Line 3317  c     $              'ETA2','ETA2',
3317                    rymm = resyPAM                    rymm = resyPAM
3318                    distmin = distance                    distmin = distance
3319                    iclm = icy                    iclm = icy
3320  c                 dedxmm = dedx(icy)  !(1)  c                 dedxmm = sgnl(icy)  !(1)
3321                    dedxmmx = 0.        !(1)                    dedxmmx = 0.        !(1)
3322                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3323                 endif                                   endif                  
3324  11882          continue  11882          continue
3325              enddo               !end loop on cluster inside couples              enddo               !end loop on cluster inside couples
3326  *----- single clusters -----------------------------------------------      *----- single clusters -----------------------------------------------  
3327    c            print*,'## ncls(',ip,') ',ncls(ip)
3328              do ic=1,ncls(ip)    !loop on single clusters              do ic=1,ncls(ip)    !loop on single clusters
3329    c               print*,'-',ic,'-'
3330                 icl=cls(ip,ic)                 icl=cls(ip,ic)
3331  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
3332                 if(cl_used(icl).ne.0.or.     !if the cluster is already used !(3)                 if(cl_used(icl).ne.0.or.     !if the cluster is already used !(3)
# Line 3512  c               if(cl_used(icl).eq.1.or. Line 3334  c               if(cl_used(icl).eq.1.or.
3334       $              .false.)goto 18882      !jump to the next singlet       $              .false.)goto 18882      !jump to the next singlet
3335                 if(mod(VIEW(icl),2).eq.0)then!<---- X view                 if(mod(VIEW(icl),2).eq.0)then!<---- X view
3336                    call xyz_PAM(icl,0,ist,                    call xyz_PAM(icl,0,ist,
 c     $                 'ETA2','ETA2',  
3337       $                 PFA,PFA,       $                 PFA,PFA,
3338       $                 AXV_STORE(nplanes-ip+1,ibest),0.)       $                 AXV_STORE(nplanes-ip+1,ibest),0.,
3339         $                 bxyz(1),
3340         $                 bxyz(2)
3341         $                 )
3342                 else                         !<---- Y view                 else                         !<---- Y view
3343                    call xyz_PAM(0,icl,ist,                    call xyz_PAM(0,icl,ist,
 c     $                 'ETA2','ETA2',  
3344       $                 PFA,PFA,       $                 PFA,PFA,
3345       $                 0.,AYV_STORE(nplanes-ip+1,ibest))       $                 0.,AYV_STORE(nplanes-ip+1,ibest),
3346         $                 bxyz(1),
3347         $                 bxyz(2)
3348         $                 )
3349                 endif                 endif
3350    
3351                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3352                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3353                 if(DEBUG)print*,'( cl-s ',icl                 if(DEBUG.EQ.1)print*,'( cl-s ',icl
3354       $              ,' ) normalized distance ',distance       $              ,' ) distance ',distance
3355                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3356    c                  if(DEBUG.EQ.1)print*,'YES'
3357                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3358                    ymm_A = yPAM_A                    ymm_A = yPAM_A
3359                    zmm_A = zPAM_A                    zmm_A = zPAM_A
# Line 3537  c     $                 'ETA2','ETA2', Line 3364  c     $                 'ETA2','ETA2',
3364                    rymm = resyPAM                    rymm = resyPAM
3365                    distmin = distance                      distmin = distance  
3366                    iclm = icl                    iclm = icl
 c                  dedxmm = dedx(icl)                   !(1)  
3367                    if(mod(VIEW(icl),2).eq.0)then !<---- X view                    if(mod(VIEW(icl),2).eq.0)then !<---- X view
3368                       dedxmmx = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmx = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3369                       dedxmmy = 0.                       !(1)                       dedxmmy = 0.                  
3370                    else          !<---- Y view                    else          !<---- Y view
3371                       dedxmmx = 0.                       !(1)                       dedxmmx = 0.                  
3372                       dedxmmy = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmy = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3373                    endif                    endif
3374                 endif                                   endif                  
3375  18882          continue  18882          continue
3376              enddo               !end loop on single clusters              enddo               !end loop on single clusters
3377    c            print*,'## distmin ', distmin,' clinc ',clinc
3378    
3379              if(distmin.le.clinc)then                    c     QUIQUI------------
3380                  c     anche qui: non ci vuole clinc???
3381                 CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<                  if(iclm.ne.0)then
 *              ----------------------------  
 c               print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'  
3382                 if(mod(VIEW(iclm),2).eq.0)then                 if(mod(VIEW(iclm),2).eq.0)then
3383                    XGOOD(nplanes-ip+1)=1.                    clincnew=
3384                    resx(nplanes-ip+1)=rxmm       $                 20*
3385                    if(DEBUG)print*,'%%%% included X-cl ',iclm       $                 sqrt(rxmm**2+RCHI2_STORE(ibest)*k(ip)*cov(1,1))
3386  c                  if(.true.)print*,'%%%% included X-cl ',iclm                 else if(mod(VIEW(iclm),2).ne.0)then
3387       $                 ,'( chi^2, ',RCHI2_STORE(ibest)                    clincnew=
3388       $                 ,', norm.dist.= ',distmin       $                 10*
3389       $                 ,', cut ',clinc,' )'       $                 sqrt(rymm**2+RCHI2_STORE(ibest)*k(ip)*cov(2,2))
3390                 else                 endif
3391                    YGOOD(nplanes-ip+1)=1.  c     QUIQUI------------
3392                    resy(nplanes-ip+1)=rymm                
3393                    if(DEBUG)print*,'%%%% included Y-cl ',iclm                 if(distmin.le.clincnew)then   !QUIQUI
3394  c                  if(.true.)print*,'%%%% included Y-cl ',iclm  c     if(distmin.le.clinc)then          !QUIQUI          
3395       $                 ,'( chi^2, ',RCHI2_STORE(ibest)                    
3396       $                 ,', norm.dist.= ', distmin                    CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<    
3397       $                 ,', cut ',clinc,' )'  *     ----------------------------
3398    c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
3399                      if(mod(VIEW(iclm),2).eq.0)then
3400                         XGOOD(nplanes-ip+1)=1.
3401                         resx(nplanes-ip+1)=rxmm
3402                         if(DEBUG.EQ.1)print*,'%%%% included X-cl ',iclm
3403         $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3404         $                    ,', dist.= ',distmin
3405         $                    ,', cut ',clinc,' )'
3406                      else
3407                         YGOOD(nplanes-ip+1)=1.
3408                         resy(nplanes-ip+1)=rymm
3409                         if(DEBUG.EQ.1)print*,'%%%% included Y-cl ',iclm
3410         $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3411         $                    ,', dist.= ', distmin
3412         $                    ,', cut ',clinc,' )'
3413                      endif
3414    c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
3415    *     ----------------------------
3416                      xm_A(nplanes-ip+1) = xmm_A
3417                      ym_A(nplanes-ip+1) = ymm_A
3418                      zm_A(nplanes-ip+1) = zmm_A
3419                      xm_B(nplanes-ip+1) = xmm_B
3420                      ym_B(nplanes-ip+1) = ymm_B
3421                      zm_B(nplanes-ip+1) = zmm_B
3422                      zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.
3423    c$$$                  zm(nplanes-ip+1) =
3424    c$$$     $                 z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
3425                      dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3426                      dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
3427    *     ----------------------------
3428                 endif                 endif
 c               print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'  
 *              ----------------------------  
                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)  
 *              ----------------------------  
3429              endif              endif
3430           endif           endif
3431   133     continue   133     continue
# Line 3600  c              dedxtrk(nplanes-ip+1) = d Line 3444  c              dedxtrk(nplanes-ip+1) = d
3444  *                                                 *  *                                                 *
3445  *                                                 *  *                                                 *
3446  **************************************************  **************************************************
 cccccc 12/08/2006 modified by elena ---> (1)  
3447  *  *
3448        subroutine clean_XYclouds(ibest,iflag)        subroutine clean_XYclouds(ibest,iflag)
3449    
3450        include 'commontracker.f'        include 'commontracker.f'
3451        include 'level1.f'        include 'level1.f'
3452        include 'common_momanhough.f'        include 'common_momanhough.f'
3453  c      include 'momanhough_init.f'        include 'level2.f'      
       include 'level2.f'        !(1)  
 c      include 'calib.f'  
 c      include 'level1.f'  
   
3454    
3455          if(DEBUG.EQ.1)print*,'clean_XYclouds:'
3456    
3457        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3458    
# Line 3622  c      include 'level1.f' Line 3462  c      include 'level1.f'
3462              if(id.ne.0)then              if(id.ne.0)then
3463                 iclx=clx(ip,icp_cp(id))                 iclx=clx(ip,icp_cp(id))
3464                 icly=cly(ip,icp_cp(id))                 icly=cly(ip,icp_cp(id))
3465  c               cl_used(iclx)=1  !tag used clusters  c$$$               cl_used(iclx)=ntrk  !tag used clusters
3466  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)  
3467              elseif(icl.ne.0)then              elseif(icl.ne.0)then
3468  c               cl_used(icl)=1   !tag used clusters  c$$$               cl_used(icl)=ntrk   !tag used clusters
                cl_used(icl)=ntrk   !tag used clusters !1)  
3469              endif              endif
3470                            
 c               if(DEBUG)then  
 c                  print*,ip,' <<< ',id  
 c               endif  
3471  *     -----------------------------  *     -----------------------------
3472  *     remove the couple from clouds  *     remove the couple from clouds
3473  *     remove also vitual couples containing the  *     remove also vitual couples containing the
# Line 3650  c               endif Line 3484  c               endif
3484       $              cly(ip,icp).eq.icl       $              cly(ip,icp).eq.icl
3485       $              )then       $              )then
3486                    id=id_cp(ip,icp,1)                    id=id_cp(ip,icp,1)
3487                    if(DEBUG)then                    if(DEBUG.EQ.1)then
3488                       print*,ip,' <<< cp ',id                       print*,ip,' <<< cp ',id
3489       $                    ,' ( cl-x '       $                    ,' ( cl-x '
3490       $                    ,clx(ip,icp)       $                    ,clx(ip,icp)
# Line 3684  c               endif Line 3518  c               endif
3518    
3519    
3520    
 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$$$  
   
3521    
3522    
3523  *     ****************************************************  *     ****************************************************
# Line 3791  c$$$ Line 3528  c$$$
3528        include 'level1.f'        include 'level1.f'
3529        include 'common_momanhough.f'        include 'common_momanhough.f'
3530        include 'level2.f'        include 'level2.f'
 c      include 'level1.f'  
3531    
3532    *     ---------------------------------
3533    *     variables initialized from level1
3534    *     ---------------------------------
3535        do i=1,nviews        do i=1,nviews
3536           good2(i)=good1(i)           good2(i)=good1(i)
3537             do j=1,nva1_view
3538                vkflag(i,j)=1
3539                if(cnnev(i,j).le.0)then
3540                   vkflag(i,j)=cnnev(i,j)
3541                endif
3542             enddo
3543        enddo        enddo
3544    *     ----------------
3545    *     level2 variables
3546    *     ----------------
3547        NTRK = 0        NTRK = 0
3548        do it=1,NTRKMAX        do it=1,NTRKMAX
3549           IMAGE(IT)=0           IMAGE(IT)=0
# Line 3808  c      include 'level1.f' Line 3554  c      include 'level1.f'
3554              ZM_nt(IP,IT) = 0              ZM_nt(IP,IT) = 0
3555              RESX_nt(IP,IT) = 0              RESX_nt(IP,IT) = 0
3556              RESY_nt(IP,IT) = 0              RESY_nt(IP,IT) = 0
3557                TAILX_nt(IP,IT) = 0
3558                TAILY_nt(IP,IT) = 0
3559                XBAD(IP,IT) = 0
3560                YBAD(IP,IT) = 0
3561              XGOOD_nt(IP,IT) = 0              XGOOD_nt(IP,IT) = 0
3562              YGOOD_nt(IP,IT) = 0              YGOOD_nt(IP,IT) = 0
3563                LS(IP,IT) = 0
3564              DEDX_X(IP,IT) = 0              DEDX_X(IP,IT) = 0
3565              DEDX_Y(IP,IT) = 0              DEDX_Y(IP,IT) = 0
3566              CLTRX(IP,IT) = 0              CLTRX(IP,IT) = 0
3567              CLTRY(IP,IT) = 0              CLTRY(IP,IT) = 0
3568                multmaxx(ip,it) = 0
3569                seedx(ip,it)    = 0
3570                xpu(ip,it)      = 0
3571                multmaxy(ip,it) = 0
3572                seedy(ip,it)    = 0
3573                ypu(ip,it)      = 0
3574           enddo           enddo
3575           do ipa=1,5           do ipa=1,5
3576              AL_nt(IPA,IT) = 0              AL_nt(IPA,IT) = 0
# Line 3833  c      include 'level1.f' Line 3590  c      include 'level1.f'
3590          ys(1,ip)=0          ys(1,ip)=0
3591          ys(2,ip)=0          ys(2,ip)=0
3592          sgnlys(ip)=0          sgnlys(ip)=0
3593            sxbad(ip)=0
3594            sybad(ip)=0
3595            multmaxsx(ip)=0
3596            multmaxsy(ip)=0
3597        enddo        enddo
3598        end        end
3599    
# Line 3865  c      include 'level1.f' Line 3626  c      include 'level1.f'
3626           alfayz1_nt(idb)=0                 alfayz1_nt(idb)=0      
3627           alfayz2_nt(idb)=0                 alfayz2_nt(idb)=0      
3628        enddo        enddo
3629        do itr=1,ntrpl_max_nt        do itr=1,ntrpt_max_nt
3630           tr_cloud_nt(itr)=0           tr_cloud_nt(itr)=0
3631           alfaxz1_nt(itr)=0                 alfaxz1_nt(itr)=0      
3632           alfaxz2_nt(itr)=0                 alfaxz2_nt(itr)=0      
# Line 3894  c      include 'level1.f' Line 3655  c      include 'level1.f'
3655          alfayz1(idb)=0                    alfayz1(idb)=0          
3656          alfayz2(idb)=0                    alfayz2(idb)=0          
3657        enddo                            enddo                    
3658        do itr=1,ntrpl_max                do itr=1,ntrpt_max        
3659          tr_cloud(itr)=0                  tr_cloud(itr)=0        
3660          cpxz1(itr)=0                      cpxz1(itr)=0            
3661          cpxz2(itr)=0                      cpxz2(itr)=0            
# Line 3942  c      include 'level1.f' Line 3703  c      include 'level1.f'
3703    
3704            
3705        include 'commontracker.f'        include 'commontracker.f'
 c      include 'level1.f'  
3706        include 'level1.f'        include 'level1.f'
3707        include 'common_momanhough.f'        include 'common_momanhough.f'
3708        include 'level2.f'        include 'level2.f'
3709        include 'common_mini_2.f'        include 'common_mini_2.f'
3710        real sinth,phi,pig              include 'calib.f'
3711    
3712          character*10 PFA
3713          common/FINALPFA/PFA
3714    
3715          real sinth,phi,pig
3716          integer ssensor,sladder
3717        pig=acos(-1.)        pig=acos(-1.)
3718    
3719    *     -------------------------------------
3720        chi2_nt(ntr)        = sngl(chi2)        chi2_nt(ntr)        = sngl(chi2)
3721        nstep_nt(ntr)       = nstep        nstep_nt(ntr)       = nstep
3722    *     -------------------------------------
3723        phi   = al(4)                  phi   = al(4)          
3724        sinth = al(3)                    sinth = al(3)            
3725        if(sinth.lt.0)then              if(sinth.lt.0)then      
# Line 3965  c      include 'level1.f' Line 3732  c      include 'level1.f'
3732       $     phi = phi + 2*pig         $     phi = phi + 2*pig  
3733        al(4) = phi                      al(4) = phi              
3734        al(3) = sinth                    al(3) = sinth            
   
3735        do i=1,5        do i=1,5
3736           al_nt(i,ntr)     = sngl(al(i))           al_nt(i,ntr)     = sngl(al(i))
3737           do j=1,5           do j=1,5
3738              coval(i,j,ntr) = sngl(cov(i,j))              coval(i,j,ntr) = sngl(cov(i,j))
3739           enddo           enddo
3740        enddo        enddo
3741          *     -------------------------------------      
3742        do ip=1,nplanes           ! loop on planes        do ip=1,nplanes           ! loop on planes
3743           xgood_nt(ip,ntr) = int(xgood(ip))           xgood_nt(ip,ntr) = int(xgood(ip))
3744           ygood_nt(ip,ntr) = int(ygood(ip))           ygood_nt(ip,ntr) = int(ygood(ip))
# Line 3981  c      include 'level1.f' Line 3747  c      include 'level1.f'
3747           zm_nt(ip,ntr)    = sngl(zm(ip))           zm_nt(ip,ntr)    = sngl(zm(ip))
3748           RESX_nt(IP,ntr)  = sngl(resx(ip))           RESX_nt(IP,ntr)  = sngl(resx(ip))
3749           RESY_nt(IP,ntr)  = sngl(resy(ip))           RESY_nt(IP,ntr)  = sngl(resy(ip))
3750             TAILX_nt(IP,ntr) = 0.
3751             TAILY_nt(IP,ntr) = 0.
3752           xv_nt(ip,ntr)    = sngl(xv(ip))           xv_nt(ip,ntr)    = sngl(xv(ip))
3753           yv_nt(ip,ntr)    = sngl(yv(ip))           yv_nt(ip,ntr)    = sngl(yv(ip))
3754           zv_nt(ip,ntr)    = sngl(zv(ip))           zv_nt(ip,ntr)    = sngl(zv(ip))
3755           axv_nt(ip,ntr)   = sngl(axv(ip))           axv_nt(ip,ntr)   = sngl(axv(ip))
3756           ayv_nt(ip,ntr)   = sngl(ayv(ip))           ayv_nt(ip,ntr)   = sngl(ayv(ip))  
3757           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)  
3758           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)             factor = sqrt(
3759         $        tan( acos(-1.) * sngl(axv(ip)) /180. )**2 +
3760         $        tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 +
3761         $        1. )
3762    
3763             dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)
3764             dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)  
3765        
3766           id  = CP_STORE(ip,IDCAND)           ax   = axv_nt(ip,ntr)
3767             ay   = ayv_nt(ip,ntr)
3768             bfx  = BX_STORE(ip,IDCAND)
3769             bfy  = BY_STORE(ip,IDCAND)
3770    c$$$         if(ip.eq.6) ax = -1. * axv_nt(ip,ntr)
3771    c$$$         if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND)
3772    c$$$         tgtemp   = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001
3773    c$$$         angx     = 180.*atan(tgtemp)/acos(-1.)
3774    c$$$         tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
3775    c$$$         angy    = 180.*atan(tgtemp)/acos(-1.)
3776    
3777             angx = effectiveangle(ax,2*ip,bfy)
3778             angy = effectiveangle(ay,2*ip-1,bfx)
3779            
3780            
3781    c         print*,'* ',ip,bfx,bfy,angx,angy
3782    
3783             id  = CP_STORE(ip,IDCAND) ! couple id
3784           icl = CLS_STORE(ip,IDCAND)           icl = CLS_STORE(ip,IDCAND)
3785             ssensor = -1
3786             sladder = -1
3787             ssensor = SENSOR_STORE(ip,IDCAND)
3788             sladder = LADDER_STORE(ip,IDCAND)
3789             if(ip.eq.6.and.ssensor.ne.0)ssensor = 3 - ssensor !notazione paolo x align
3790             LS(IP,ntr)      = ssensor+10*sladder
3791    
3792           if(id.ne.0)then           if(id.ne.0)then
3793    c           >>> is a couple
3794              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))
3795              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))
3796  c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)  
3797                cl_used(cltrx(ip,ntr)) = 1     !tag used clusters          
3798                cl_used(cltry(ip,ntr)) = 1     !tag used clusters          
3799    
3800                xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))
3801                ybad(ip,ntr)= nbadstrips(4,cly(nplanes-ip+1,icp_cp(id)))
3802    
3803    
3804                if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0)
3805         $           dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3806                if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)
3807         $           dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3808    
3809                multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
3810         $                         +10000*mult(cltrx(ip,ntr))
3811                seedx(ip,ntr)    = clsignal(indmax(cltrx(ip,ntr)))
3812         $           /clsigma(indmax(cltrx(ip,ntr)))
3813                call applypfa(PFA,cltrx(ip,ntr),angx,corr,res)
3814                xpu(ip,ntr)      = corr
3815    
3816                multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
3817         $                         +10000*mult(cltry(ip,ntr))
3818                seedy(ip,ntr)    = clsignal(indmax(cltry(ip,ntr)))
3819         $           /clsigma(indmax(cltry(ip,ntr)))
3820                call applypfa(PFA,cltry(ip,ntr),angy,corr,res)
3821                ypu(ip,ntr)      = corr
3822    
3823           elseif(icl.ne.0)then           elseif(icl.ne.0)then
3824              if(mod(VIEW(icl),2).eq.0)cltrx(ip,ntr)=icl  
3825              if(mod(VIEW(icl),2).eq.1)cltry(ip,ntr)=icl  
3826  c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)              cl_used(icl) = 1    !tag used clusters          
3827    
3828                if(mod(VIEW(icl),2).eq.0)then
3829                   cltrx(ip,ntr)=icl
3830                   xbad(ip,ntr) = nbadstrips(4,icl)
3831    
3832                   if(nsatstrips(icl).gt.0)dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3833    
3834                   multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
3835         $                         +10000*mult(cltrx(ip,ntr))
3836                   seedx(ip,ntr)    = clsignal(indmax(cltrx(ip,ntr)))
3837         $           /clsigma(indmax(cltrx(ip,ntr)))
3838                   call applypfa(PFA,cltrx(ip,ntr),angx,corr,res)
3839                   xpu(ip,ntr)      = corr
3840    
3841                elseif(mod(VIEW(icl),2).eq.1)then
3842                   cltry(ip,ntr)=icl
3843                   ybad(ip,ntr) = nbadstrips(4,icl)
3844    
3845                   if(nsatstrips(icl).gt.0)dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3846    
3847                   multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
3848         $                         +10000*mult(cltry(ip,ntr))
3849                   seedy(ip,ntr)    = clsignal(indmax(cltry(ip,ntr)))
3850         $           /clsigma(indmax(cltry(ip,ntr)))
3851                   call applypfa(PFA,cltry(ip,ntr),angy,corr,res)
3852                   ypu(ip,ntr)      = corr
3853                  
3854                endif
3855    
3856           endif                     endif          
3857    
3858        enddo        enddo
3859    
3860          if(DEBUG.eq.1)then
3861             print*,'> STORING TRACK ',ntr
3862             print*,'clusters: '
3863             do ip=1,6
3864                print*,'> ',ip,' -- ',cltrx(ip,ntr),cltry(ip,ntr)
3865             enddo
3866          endif
3867    
3868    c$$$      print*,(xgood(i),i=1,6)
3869    c$$$      print*,(ygood(i),i=1,6)
3870    c$$$      print*,(ls(i,ntr),i=1,6)
3871    c$$$      print*,(dedx_x(i,ntr),i=1,6)
3872    c$$$      print*,(dedx_y(i,ntr),i=1,6)
3873    c$$$      print*,'-----------------------'
3874    
3875        end        end
3876    
# Line 4015  c            print*,ip,' ',cltrx(ip,ntr) Line 3883  c            print*,ip,' ',cltrx(ip,ntr)
3883  *     -------------------------------------------------------  *     -------------------------------------------------------
3884    
3885        include 'commontracker.f'        include 'commontracker.f'
 c      include 'level1.f'  
3886        include 'calib.f'        include 'calib.f'
3887        include 'level1.f'        include 'level1.f'
3888        include 'common_momanhough.f'        include 'common_momanhough.f'
# Line 4023  c      include 'level1.f' Line 3890  c      include 'level1.f'
3890        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3891    
3892  *     count #cluster per plane not associated to any track  *     count #cluster per plane not associated to any track
 c      good2=1!.true.  
3893        nclsx = 0        nclsx = 0
3894        nclsy = 0        nclsy = 0
3895    
3896        do iv = 1,nviews        do iv = 1,nviews
3897           if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv)  c         if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv)
3898             good2(iv) = good2(iv) + mask_view(iv)*2**8
3899        enddo        enddo
3900    
3901          if(DEBUG.eq.1)then
3902             print*,'> STORING SINGLETS '
3903          endif
3904    
3905        do icl=1,nclstr1        do icl=1,nclstr1
3906    
3907             ip=nplanes-npl(VIEW(icl))+1            
3908            
3909           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
3910              ip=nplanes-npl(VIEW(icl))+1              
3911    cc            print*,' ic ',icl,' --- stored '
3912    
3913              if(mod(VIEW(icl),2).eq.0)then !=== X views              if(mod(VIEW(icl),2).eq.0)then !=== X views
3914    
3915                 nclsx = nclsx + 1                 nclsx = nclsx + 1
3916                 planex(nclsx) = ip                 planex(nclsx) = ip
3917                 sgnlxs(nclsx) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3918                   if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx)
3919                 clsx(nclsx)   = icl                 clsx(nclsx)   = icl
3920                   sxbad(nclsx)  = nbadstrips(1,icl)
3921                   multmaxsx(nclsx) = maxs(icl)+10000*mult(icl)
3922                  
3923    cc               print*,icl,' >>>> ',sxbad(nclsx)
3924    
3925                 do is=1,2                 do is=1,2
3926  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
3927                    call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)  c                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
3928                      call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.,0.,0.)
3929                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2
3930                 enddo                 enddo
3931  c$$$               print*,'nclsx         ',nclsx  c$$$               print*,'nclsx         ',nclsx
# Line 4052  c$$$               print*,'xs(2,nclsx)   Line 3936  c$$$               print*,'xs(2,nclsx)  
3936              else                          !=== Y views              else                          !=== Y views
3937                 nclsy = nclsy + 1                 nclsy = nclsy + 1
3938                 planey(nclsy) = ip                 planey(nclsy) = ip
3939                 sgnlys(nclsy) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3940                   if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy)
3941                 clsy(nclsy)   = icl                 clsy(nclsy)   = icl
3942                   sybad(nclsy)  = nbadstrips(1,icl)
3943                   multmaxsy(nclsy) = maxs(icl)+10000*mult(icl)
3944    
3945    cc               print*,icl,' >>>> ',sybad(nclsy)
3946    
3947                 do is=1,2                 do is=1,2
3948  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
3949                    call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)  c                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
3950                      call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.,0.,0.)
3951                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2
3952                 enddo                 enddo
3953  c$$$               print*,'nclsy         ',nclsy  c$$$               print*,'nclsy         ',nclsy
# Line 4066  c$$$               print*,'ys(1,nclsy)   Line 3957  c$$$               print*,'ys(1,nclsy)  
3957  c$$$               print*,'ys(2,nclsy)   ',ys(2,nclsy)  c$$$               print*,'ys(2,nclsy)   ',ys(2,nclsy)
3958              endif              endif
3959           endif           endif
 c      print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip)  
3960    
3961  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO
3962           whichtrack(icl) = cl_used(icl)           whichtrack(icl) = cl_used(icl)
3963    *     --------------------------------------------------
3964    *     per non perdere la testa...
3965    *     whichtrack(icl) e` una variabile del common level1
3966    *     che serve solo per sapere quali cluster sono stati
3967    *     associati ad una traccia, e permettere di salvare
3968    *     solo questi nell'albero di uscita
3969    *     --------------------------------------------------
3970                    
3971        enddo        enddo
3972        end        end
3973    
Legend:
Removed from v.1.11  
changed lines
  Added in v.1.36
  ViewVC Help
Powered by ViewVC 1.1.23