/[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.16 by pam-fi, Thu Nov 30 17:04:27 2006 UTC revision 1.32 by bongi, Tue Sep 4 15:44:49 2007 UTC
# Line 20  Line 20 
20        include 'calib.f'        include 'calib.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 41  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 75  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 182  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 239  c$$$         if(ibest.eq.0)goto 880 !>> Line 273  c$$$         if(ibest.eq.0)goto 880 !>>
273  *     2nd) increasing chi**2  *     2nd) increasing chi**2
274  *     -------------------------------------------------------  *     -------------------------------------------------------
275           rchi2best=1000000000.           rchi2best=1000000000.
276           ndofbest=0             !(1)           ndofbest=0            
277           do i=1,ntracks           do i=1,ntracks
278             ndof=0               !(1)             ndof=0              
279             do ii=1,nplanes      !(1)             do ii=1,nplanes    
280               ndof=ndof          !(1)               ndof=ndof        
281       $            +int(xgood_store(ii,i)) !(1)       $            +int(xgood_store(ii,i))
282       $            +int(ygood_store(ii,i)) !(1)       $            +int(ygood_store(ii,i))
283             enddo                !(1)             enddo              
284             if(ndof.gt.ndofbest)then !(1)             if(ndof.gt.ndofbest)then
285               ibest=i               ibest=i
286               rchi2best=RCHI2_STORE(i)               rchi2best=RCHI2_STORE(i)
287               ndofbest=ndof      !(1)               ndofbest=ndof    
288             elseif(ndof.eq.ndofbest)then !(1)             elseif(ndof.eq.ndofbest)then
289               if(RCHI2_STORE(i).lt.rchi2best.and.               if(RCHI2_STORE(i).lt.rchi2best.and.
290       $            RCHI2_STORE(i).gt.0)then       $            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    
# Line 298  c$$$         enddo Line 332  c$$$         enddo
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 314  c$$$         enddo Line 350  c$$$         enddo
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 324  c$$$         enddo Line 361  c$$$         enddo
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(VERBOSE)iprint=1           if(VERBOSE.EQ.1)iprint=1
366           if(DEBUG)iprint=2           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(VERBOSE)then              if(VERBOSE.EQ.1)then
370                 print *,                 print *,
371       $              '*** MINIMIZATION FAILURE *** (after refinement) '       $              '*** MINIMIZATION FAILURE *** (after refinement) '
372       $              ,iev       $              ,iev
# Line 344  c$$$               print*,'------------- Line 381  c$$$               print*,'-------------
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 355  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 371  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 402  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 438  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 476  c     $        rchi2best.lt.15..and. Line 628  c     $        rchi2best.lt.15..and.
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 514  c     $        rchi2best.lt.15..and. Line 668  c     $        rchi2best.lt.15..and.
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    c      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    c$$$      print*,'## stripx,stripy ',stripx,stripy
809    
         
810  c===========================================================  c===========================================================
811  C     COUPLE  C     COUPLE
812  C===========================================================  C===========================================================
# Line 697  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 738  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 764  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 790  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 813  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 865  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 881  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         $           ,' does not exists (nclstr1=',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             xm(ip) = -100.
1131             ym(ip) = -100.
1132             zm(ip) = (zPAM_A+zPAM_B)/2.
1133             xm_A(ip) = xPAM_A
1134             ym_A(ip) = yPAM_A
1135             xm_B(ip) = xPAM_B
1136             ym_B(ip) = yPAM_B
1137    
1138    c         zv(ip) = (zPAM_A+zPAM_B)/2.
1139            
1140          elseif(icx.ne.0.and.icy.eq.0)then
1141    
1142             ipx=npl(VIEW(icx))
1143    c$$$         if((nplanes-ipx+1).ne.ip)
1144    c$$$     $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1145    c$$$     $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1146    
1147             if( (nplanes-ipx+1).ne.ip )then
1148                print*,'xyzpam: ***WARNING*** cluster ',icx
1149         $           ,' does not belong to plane: ',ip
1150                icx = -1*icx
1151                return
1152             endif
1153    
1154             call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1155          
1156             xgood(ip) = 1.
1157             ygood(ip) = 0.
1158             resx(ip)  = resxPAM
1159             resy(ip)  = 1000.
1160    
1161             xm(ip) = -100.
1162             ym(ip) = -100.
1163             zm(ip) = (zPAM_A+zPAM_B)/2.
1164             xm_A(ip) = xPAM_A
1165             ym_A(ip) = yPAM_A
1166             xm_B(ip) = xPAM_B
1167             ym_B(ip) = yPAM_B
1168            
1169    c         zv(ip) = (zPAM_A+zPAM_B)/2.
1170    
1171          else
1172    
1173             il = 2
1174             if(lad.ne.0)il=lad
1175             is = 1
1176             if(sensor.ne.0)is=sensor
1177    c         print*,nplanes-ip+1,il,is
1178    
1179             xgood(ip) = 0.
1180             ygood(ip) = 0.
1181             resx(ip)  = 1000.
1182             resy(ip)  = 1000.
1183    
1184             xm(ip) = -100.
1185             ym(ip) = -100.          
1186             zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1187             xm_A(ip) = 0.
1188             ym_A(ip) = 0.
1189             xm_B(ip) = 0.
1190             ym_B(ip) = 0.
1191    
1192    c         zv(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1193    
1194          endif
1195    
1196          if(DEBUG.EQ.1)then
1197    c         print*,'----------------------------- track coord'
1198    22222    format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
1199             write(*,22222)ip,zm(ip),xm(ip),ym(ip)
1200         $        ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip)
1201         $        ,xgood(ip),ygood(ip),resx(ip),resy(ip)
1202    c$$$         print*,'-----------------------------'
1203          endif
1204          end
1205    
1206  ********************************************************************************  ********************************************************************************
1207  ********************************************************************************  ********************************************************************************
# Line 912  c         print*,'A-(',xPAM_A,yPAM_A,') Line 1223  c         print*,'A-(',xPAM_A,yPAM_A,')
1223  *      *    
1224  ********************************************************************************  ********************************************************************************
1225    
1226        real function distance_to(XPP,YPP)        real function distance_to(rXPP,rYPP)
1227    
1228        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
1229    
# Line 921  c         print*,'A-(',xPAM_A,yPAM_A,') Line 1232  c         print*,'A-(',xPAM_A,yPAM_A,')
1232  *     ( i.e. distance/resolution )  *     ( i.e. distance/resolution )
1233  *     -----------------------------------  *     -----------------------------------
1234    
1235          real rXPP,rYPP
1236          double precision XPP,YPP
1237        double precision distance,RE        double precision distance,RE
1238        double precision BETA,ALFA,xmi,ymi        double precision BETA,ALFA,xmi,ymi
1239    
1240          XPP=DBLE(rXPP)
1241          YPP=DBLE(rYPP)
1242    
1243  *     ----------------------  *     ----------------------
1244        if (        if (
1245       +     xPAM.eq.0.and.       +     xPAM.eq.0.and.
# Line 966  c         print*,'A-(',xPAM_A,yPAM_A,') Line 1282  c         print*,'A-(',xPAM_A,yPAM_A,')
1282           endif                   endif        
1283    
1284           distance=           distance=
1285       $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2       $       ((xmi-XPP)**2+(ymi-YPP)**2)!QUIQUI
1286    cc     $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2
1287           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1288    
1289  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 991  c$$$         print*,' resolution ',re Line 1308  c$$$         print*,' resolution ',re
1308  *     ----------------------  *     ----------------------
1309                    
1310           distance=           distance=
1311       $        ((xPAM-XPP)/resxPAM)**2       $       ((xPAM-XPP))**2 !QUIQUI
1312       $        +       $       +
1313       $        ((yPAM-YPP)/resyPAM)**2       $       ((yPAM-YPP))**2
1314    c$$$     $        ((xPAM-XPP)/resxPAM)**2
1315    c$$$     $        +
1316    c$$$     $        ((yPAM-YPP)/resyPAM)**2
1317           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1318    
1319  c$$$         print*,xPAM,yPAM,zPAM  c$$$         print*,xPAM,yPAM,zPAM
# Line 1002  c$$$         print*,' resolution ',resxP Line 1322  c$$$         print*,' resolution ',resxP
1322                    
1323        else        else
1324                    
1325           print*  c         print*
1326       $        ,' function distance_to ---> wrong usage!!!'  c     $        ,' function distance_to ---> wrong usage!!!'
1327           print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM  c         print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM
1328           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 '
1329       $        ,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
1330        endif          endif  
1331    
1332        distance_to = sngl(distance)        distance_to = sngl(distance)
# Line 1046  c$$$         print*,' resolution ',resxP Line 1366  c$$$         print*,' resolution ',resxP
1366        data c1/1.,0.,0.,1./        data c1/1.,0.,0.,1./
1367        data c2/1.,-1.,-1.,1./        data c2/1.,-1.,-1.,1./
1368        data c3/1.,1.,0.,0./        data c3/1.,1.,0.,0./
1369        real*8 yvvv,xvvv        double precision yvvv,xvvv
1370        double precision xi,yi,zi        double precision xi,yi,zi
1371        double precision xrt,yrt,zrt        double precision xrt,yrt,zrt
1372        real AA,BB        real AA,BB
1373        real yvv(4),xvv(4)        double precision yvv(4),xvv(4)
1374    
1375  *     tollerance to consider the track inside the sensitive area  *     tollerance to consider the track inside the sensitive area
1376        real ptoll        real ptoll
1377        data ptoll/150./          !um        data ptoll/150./          !um
1378    
1379        external nviewx,nviewy,acoordsi,dcoord        external nviewx,nviewy,dcoordsi,dcoord
1380    
1381        nplpt = nplPAM            !plane        nplpt = nplPAM            !plane
1382        viewx = nviewx(nplpt)        viewx = nviewx(nplpt)
# Line 1071  c$$$         print*,' resolution ',resxP Line 1391  c$$$         print*,' resolution ',resxP
1391  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1392  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
1393  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1394                 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)  c               if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
1395       $              .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...
1396  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...
1397                    print*,'whichsensor: ',  c                  print*,'whichsensor: ',
1398       $                ' WARNING: false X strip: strip ',stripx  c     $                ' WARNING: false X strip: strip ',stripx
1399                 endif  c               endif
1400                 xi = acoordsi(stripx,viewx)                 xi = dcoordsi(stripx,viewx)
1401                 yi = acoordsi(stripy,viewy)                 yi = dcoordsi(stripy,viewy)
1402                 zi = 0.                 zi = 0.D0
1403  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1404  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
1405  c------------------------------------------------------------------------  c------------------------------------------------------------------------
# Line 1230  c      include 'common_analysis.f' Line 1550  c      include 'common_analysis.f'
1550        is_cp=0        is_cp=0
1551        if(id.lt.0)is_cp=1        if(id.lt.0)is_cp=1
1552        if(id.gt.0)is_cp=2        if(id.gt.0)is_cp=2
1553        if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'  c      if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'
1554    
1555        return        return
1556        end        end
# Line 1329  c      include 'level1.f' Line 1649  c      include 'level1.f'
1649        integer iflag        integer iflag
1650    
1651        integer badseed,badclx,badcly        integer badseed,badclx,badcly
1652        
1653          if(DEBUG.EQ.1)print*,'cl_to_couples:'
1654    
1655  *     init variables  *     init variables
1656        ncp_tot=0        ncp_tot=0
# Line 1359  c      include 'level1.f' Line 1681  c      include 'level1.f'
1681  *     mask views with too many clusters  *     mask views with too many clusters
1682        do iv=1,nviews        do iv=1,nviews
1683           if( ncl_view(iv).gt. nclusterlimit)then           if( ncl_view(iv).gt. nclusterlimit)then
1684              mask_view(iv) = 1  c            mask_view(iv) = 1
1685              if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '              mask_view(iv) = mask_view(iv) + 2**0
1686       $           ,nclusterlimit,' on view ', iv,' --> masked!'              if(DEBUG.EQ.1)
1687         $        print*,' * WARNING * cl_to_couple: n.clusters > '
1688         $        ,nclusterlimit,' on view ', iv,' --> masked!'
1689           endif           endif
1690        enddo        enddo
1691    
# Line 1378  c      include 'level1.f' Line 1702  c      include 'level1.f'
1702  *     ----------------------------------------------------  *     ----------------------------------------------------
1703  *     cut on charge (X VIEW)  *     cut on charge (X VIEW)
1704  *     ----------------------------------------------------  *     ----------------------------------------------------
1705           if(dedx(icx).lt.dedx_x_min)then           if(sgnl(icx).lt.dedx_x_min)then
1706              cl_single(icx)=0              cl_single(icx)=0
1707              goto 10              goto 10
1708           endif           endif
# Line 1428  c     endif Line 1752  c     endif
1752  *     ----------------------------------------------------  *     ----------------------------------------------------
1753  *     cut on charge (Y VIEW)  *     cut on charge (Y VIEW)
1754  *     ----------------------------------------------------  *     ----------------------------------------------------
1755              if(dedx(icy).lt.dedx_y_min)then              if(sgnl(icy).lt.dedx_y_min)then
1756                 cl_single(icy)=0                 cl_single(icy)=0
1757                 goto 20                 goto 20
1758              endif              endif
# Line 1474  c     endif Line 1798  c     endif
1798  *     charge correlation  *     charge correlation
1799  *     (modified to be applied only below saturation... obviously)  *     (modified to be applied only below saturation... obviously)
1800    
1801                 if(  .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx)                 if(  .not.(sgnl(icy).gt.chsaty.and.sgnl(icx).gt.chsatx)
1802       $              .and.       $              .and.
1803       $              .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx)       $              .not.(sgnl(icy).lt.chmipy.and.sgnl(icx).lt.chmipx)
1804       $              .and.       $              .and.
1805       $              (badclx.eq.1.and.badcly.eq.1)       $              (badclx.eq.1.and.badcly.eq.1)
1806       $              .and.       $              .and.
1807       $              .true.)then       $              .true.)then
1808    
1809                    ddd=(dedx(icy)                    ddd=(sgnl(icy)
1810       $                 -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))       $                 -kch(nplx,nldx)*sgnl(icx)-cch(nplx,nldx))
1811                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
1812    
1813  c                  cut = chcut * sch(nplx,nldx)  c                  cut = chcut * sch(nplx,nldx)
1814    
1815                    sss=(kch(nplx,nldx)*dedx(icy)+dedx(icx)                    sss=(kch(nplx,nldx)*sgnl(icy)+sgnl(icx)
1816       $                 -kch(nplx,nldx)*cch(nplx,nldx))       $                 -kch(nplx,nldx)*cch(nplx,nldx))
1817                    sss=sss/sqrt(kch(nplx,nldx)**2+1)                    sss=sss/sqrt(kch(nplx,nldx)**2+1)
1818                    cut = chcut * (16 + sss/50.)                    cut = chcut * (16 + sss/50.)
# Line 1499  c                  cut = chcut * sch(npl Line 1823  c                  cut = chcut * sch(npl
1823                 endif                 endif
1824    
1825                 if(ncp_plane(nplx).gt.ncouplemax)then                 if(ncp_plane(nplx).gt.ncouplemax)then
1826                    if(verbose)print*,                    if(verbose.eq.1)print*,
1827       $                 '** warning ** number of identified '//       $                 '** warning ** number of identified '//
1828       $                 'couples on plane ',nplx,       $                 'couples on plane ',nplx,
1829       $                 'exceeds vector dimention '       $                 'exceeds vector dimention '
1830       $                 ,'( ',ncouplemax,' ) --> masked!'       $                 ,'( ',ncouplemax,' ) --> masked!'
1831                    mask_view(nviewx(nplx)) = 2  c                  mask_view(nviewx(nplx)) = 2
1832                    mask_view(nviewy(nply)) = 2  c                  mask_view(nviewy(nply)) = 2
1833                      mask_view(nviewx(nplx))= mask_view(nviewx(nplx))+ 2**1
1834                      mask_view(nviewy(nply))= mask_view(nviewy(nply))+ 2**1
1835                    goto 10                    goto 10
1836                 endif                 endif
1837                                
# Line 1535  c                  cut = chcut * sch(npl Line 1861  c                  cut = chcut * sch(npl
1861        enddo        enddo
1862                
1863                
1864        if(DEBUG)then        if(DEBUG.EQ.1)then
1865           print*,'clusters  ',nclstr1           print*,'clusters  ',nclstr1
1866           print*,'good    ',(cl_good(i),i=1,nclstr1)           print*,'good    ',(cl_good(i),i=1,nclstr1)
1867           print*,'singles ',(cl_single(i),i=1,nclstr1)           print*,'singlets',(cl_single(i),i=1,nclstr1)
1868           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
1869        endif        endif
1870                
# Line 1598  c      double precision xm3,ym3,zm3 Line 1924  c      double precision xm3,ym3,zm3
1924        real xc,zc,radius        real xc,zc,radius
1925  *     -----------------------------  *     -----------------------------
1926    
1927          if(DEBUG.EQ.1)print*,'cp_to_doubtrip:'
1928    
1929  *     --------------------------------------------  *     --------------------------------------------
1930  *     put a limit to the maximum number of couples  *     put a limit to the maximum number of couples
# Line 1606  c      double precision xm3,ym3,zm3 Line 1933  c      double precision xm3,ym3,zm3
1933  *     --------------------------------------------  *     --------------------------------------------
1934        do ip=1,nplanes        do ip=1,nplanes
1935           if(ncp_plane(ip).gt.ncouplelimit)then           if(ncp_plane(ip).gt.ncouplelimit)then
1936              mask_view(nviewx(ip)) = 8  c            mask_view(nviewx(ip)) = 8
1937              mask_view(nviewy(ip)) = 8  c            mask_view(nviewy(ip)) = 8
1938                mask_view(nviewx(ip)) = mask_view(nviewx(ip)) + 2**7
1939                mask_view(nviewy(ip)) = mask_view(nviewy(ip)) + 2**7
1940           endif           endif
1941        enddo        enddo
1942    
# Line 1623  c      double precision xm3,ym3,zm3 Line 1952  c      double precision xm3,ym3,zm3
1952                 icx1=clx(ip1,icp1)                 icx1=clx(ip1,icp1)
1953                 icy1=cly(ip1,icp1)                 icy1=cly(ip1,icp1)
1954  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)
1955                 call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)  c               call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)
1956                   call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.)
1957                 xm1=xPAM                 xm1=xPAM
1958                 ym1=yPAM                 ym1=yPAM
1959                 zm1=zPAM                                   zm1=zPAM                  
# Line 1632  c     print*,'***',is1,xm1,ym1,zm1 Line 1962  c     print*,'***',is1,xm1,ym1,zm1
1962                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
1963                    if(  mask_view(nviewx(ip2)).ne.0 .or.                    if(  mask_view(nviewx(ip2)).ne.0 .or.
1964       $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane       $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane
1965                    do is2=1,2    !loop on sensors -ndblt COPPIA 2                    do is2=1,2    !loop on sensors -ndblt COPPIA 2                    
                       
1966                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
1967                          icx2=clx(ip2,icp2)                          icx2=clx(ip2,icp2)
1968                          icy2=cly(ip2,icp2)                          icy2=cly(ip2,icp2)
1969  c                        call xyz_PAM  c                        call xyz_PAM
1970  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)
1971    c                        call xyz_PAM
1972    c     $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)
1973                          call xyz_PAM                          call xyz_PAM
1974       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.)
1975                          xm2=xPAM                          xm2=xPAM
1976                          ym2=yPAM                          ym2=yPAM
1977                          zm2=zPAM                          zm2=zPAM
# Line 1650  c     $                       (icx2,icy2 Line 1981  c     $                       (icx2,icy2
1981  *     (2 couples needed)  *     (2 couples needed)
1982  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1983                          if(ndblt.eq.ndblt_max)then                          if(ndblt.eq.ndblt_max)then
1984                             if(verbose)print*,                             if(verbose.eq.1)print*,
1985       $                          '** warning ** number of identified '//       $                          '** warning ** number of identified '//
1986       $                          'doublets exceeds vector dimention '       $                          'doublets exceeds vector dimention '
1987       $                          ,'( ',ndblt_max,' )'       $                          ,'( ',ndblt_max,' )'
1988  c                           good2=.false.  c                           good2=.false.
1989  c                           goto 880 !fill ntp and go to next event  c                           goto 880 !fill ntp and go to next event
1990                             do iv=1,12                             do iv=1,12
1991                                mask_view(iv) = 3  c                              mask_view(iv) = 3
1992                                  mask_view(iv) = mask_view(iv)+ 2**2
1993                             enddo                             enddo
1994                             iflag=1                             iflag=1
1995                             return                             return
# Line 1704  c$$$ Line 2036  c$$$
2036                                   icy3=cly(ip3,icp3)                                   icy3=cly(ip3,icp3)
2037  c                                 call xyz_PAM  c                                 call xyz_PAM
2038  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)
2039    c                                 call xyz_PAM
2040    c     $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)
2041                                   call xyz_PAM                                   call xyz_PAM
2042       $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)       $                                (icx3,icy3,is3,PFAdef,PFAdef
2043         $                                ,0.,0.,0.,0.)
2044                                   xm3=xPAM                                   xm3=xPAM
2045                                   ym3=yPAM                                   ym3=yPAM
2046                                   zm3=zPAM                                   zm3=zPAM
2047  *     find the circle passing through the three points  *     find the circle passing through the three points
2048    c$$$                                 call tricircle(3,xp,zp,angp,resp,chi
2049    c$$$     $                                ,xc,zc,radius,iflag)
2050                                     iflag = DEBUG
2051                                   call tricircle(3,xp,zp,angp,resp,chi                                   call tricircle(3,xp,zp,angp,resp,chi
2052       $                                ,xc,zc,radius,iflag)       $                                ,xc,zc,radius,iflag)
2053  c     print*,xc,zc,radius  c     print*,xc,zc,radius
# Line 1726  c     $                                 Line 2064  c     $                                
2064  *     (3 couples needed)  *     (3 couples needed)
2065  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2066                                   if(ntrpt.eq.ntrpt_max)then                                   if(ntrpt.eq.ntrpt_max)then
2067                                      if(verbose)print*,                                      if(verbose.eq.1)print*,
2068       $                     '** warning ** number of identified '//       $                     '** warning ** number of identified '//
2069       $                     'triplets exceeds vector dimention '       $                     'triplets exceeds vector dimention '
2070       $                    ,'( ',ntrpt_max,' )'       $                    ,'( ',ntrpt_max,' )'
2071  c                                    good2=.false.  c                                    good2=.false.
2072  c                                    goto 880 !fill ntp and go to next event  c                                    goto 880 !fill ntp and go to next event
2073                                      do iv=1,nviews                                      do iv=1,nviews
2074                                         mask_view(iv) = 4  c                                       mask_view(iv) = 4
2075                                           mask_view(iv)=mask_view(iv)+ 2**3
2076                                      enddo                                      enddo
2077                                      iflag=1                                      iflag=1
2078                                      return                                      return
# Line 1787  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp Line 2126  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp
2126   10   continue   10   continue
2127        enddo                     !end loop on planes  - COPPIA 1        enddo                     !end loop on planes  - COPPIA 1
2128                
2129        if(DEBUG)then        if(DEBUG.EQ.1)then
2130           print*,'--- doublets ',ndblt           print*,'--- doublets ',ndblt
2131           print*,'--- triplets ',ntrpt           print*,'--- triplets ',ntrpt
2132        endif        endif
# Line 1834  c      include 'momanhough_init.f' Line 2173  c      include 'momanhough_init.f'
2173        integer cp_useds1(ncouplemaxtot) ! sensor 1        integer cp_useds1(ncouplemaxtot) ! sensor 1
2174        integer cp_useds2(ncouplemaxtot) ! sensor 2        integer cp_useds2(ncouplemaxtot) ! sensor 2
2175    
2176          if(DEBUG.EQ.1)print*,'doub_to_YZcloud:'
2177    
2178  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2179  *     classification of DOUBLETS  *     classification of DOUBLETS
# Line 1960  c         if(ncpused.lt.ncpyz_min)goto 2 Line 2300  c         if(ncpused.lt.ncpyz_min)goto 2
2300  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2301    
2302           if(nclouds_yz.ge.ncloyz_max)then           if(nclouds_yz.ge.ncloyz_max)then
2303              if(verbose)print*,              if(verbose.eq.1)print*,
2304       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2305       $           'YZ clouds exceeds vector dimention '       $           'YZ clouds exceeds vector dimention '
2306       $           ,'( ',ncloyz_max,' )'       $           ,'( ',ncloyz_max,' )'
2307  c               good2=.false.  c               good2=.false.
2308  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2309              do iv=1,nviews              do iv=1,nviews
2310                 mask_view(iv) = 5  c               mask_view(iv) = 5
2311                   mask_view(iv) = mask_view(iv) + 2**4
2312              enddo              enddo
2313              iflag=1              iflag=1
2314              return              return
# Line 1987  c     ptcloud_yz_nt(nclouds_yz)=npt Line 2328  c     ptcloud_yz_nt(nclouds_yz)=npt
2328  c     print*,'>> ',ipt,db_all(ipt)  c     print*,'>> ',ipt,db_all(ipt)
2329           enddo             enddo  
2330           npt_tot=npt_tot+npt           npt_tot=npt_tot+npt
2331           if(DEBUG)then           if(DEBUG.EQ.1)then
2332              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
2333              print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'              print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'
2334              print*,'- alfayz1 ',alfayz1_av(nclouds_yz)              print*,'- alfayz1  ',alfayz1_av(nclouds_yz)
2335              print*,'- alfayz2 ',alfayz2_av(nclouds_yz)              print*,'- alfayz2  ',alfayz2_av(nclouds_yz)
2336              print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)              print*,'cp_useds1  ',(cp_useds1(icp),icp=1,ncp_tot)
2337              print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)              print*,'cp_useds2  ',(cp_useds2(icp),icp=1,ncp_tot)
2338              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)              print*,'cpcloud_yz '
2339         $           ,(cpcloud_yz(nclouds_yz,icp),icp=1,ncp_tot)
2340                print*,'hit_plane  ',(hit_plane(ip),ip=1,nplanes)
2341  c$$$            print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = '  c$$$            print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = '
2342  c$$$     $           ,ptcloud_yz(nclouds_yz)  c$$$     $           ,ptcloud_yz(nclouds_yz)
2343  c$$$            print*,'nt-uple: db_cloud(...) = '  c$$$            print*,'nt-uple: db_cloud(...) = '
# Line 2012  c$$$     $           ,(db_cloud(iii),iii Line 2355  c$$$     $           ,(db_cloud(iii),iii
2355          goto 90                          goto 90                
2356        endif                            endif                    
2357                
2358        if(DEBUG)then        if(DEBUG.EQ.1)then
2359           print*,'---------------------- '           print*,'---------------------- '
2360           print*,'Y-Z total clouds ',nclouds_yz           print*,'Y-Z total clouds ',nclouds_yz
2361           print*,' '           print*,' '
# Line 2061  c      include 'momanhough_init.f' Line 2404  c      include 'momanhough_init.f'
2404        integer cp_useds1(ncouplemaxtot) ! sensor 1        integer cp_useds1(ncouplemaxtot) ! sensor 1
2405        integer cp_useds2(ncouplemaxtot) ! sensor 2        integer cp_useds2(ncouplemaxtot) ! sensor 2
2406    
2407          if(DEBUG.EQ.1)print*,'trip_to_XZcloud:'
2408    
2409  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2410  *     classification of TRIPLETS  *     classification of TRIPLETS
2411  *     according to distance in parameter space  *     according to distance in parameter space
# Line 2118  c         tr_temp(1)=itr1 Line 2463  c         tr_temp(1)=itr1
2463       $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2                     $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2              
2464                 distance = sqrt(distance)                 distance = sqrt(distance)
2465                                
2466                 if(distance.lt.cutdistxz)then  *     ------------------------------------------------------------------------
2467    *     FORCE INCLUSION OF TRIPLETS COMPOSED BY SAME COUPLES, IGNORING THE IMAGE
2468    *     ------------------------------------------------------------------------
2469    *     (added in august 2007)
2470                   istrimage=0
2471                   if(
2472         $              abs(cpxz1(itrref)).eq.abs(cpxz1(itr2)).and.
2473         $              abs(cpxz2(itrref)).eq.abs(cpxz2(itr2)).and.
2474         $              abs(cpxz3(itrref)).eq.abs(cpxz3(itr2)).and.
2475         $              .true.)istrimage=1
2476    
2477                   if(distance.lt.cutdistxz.or.istrimage.eq.1)then
2478  c     print*,idb1,idb2,distance,' cloud ',nclouds_yz  c     print*,idb1,idb2,distance,' cloud ',nclouds_yz
2479                    if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1                    if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1
2480                    if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1                    if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1
# Line 2177  c     print*,'check cp_used' Line 2533  c     print*,'check cp_used'
2533           enddo           enddo
2534  c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet  c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet
2535           if(npt.lt.nptxz_min)goto 22288     !next triplet           if(npt.lt.nptxz_min)goto 22288     !next triplet
2536           if(nplused.lt.nplxz_min)goto 22288 !next doublet           if(nplused.lt.nplxz_min)goto 22288 !next triplet
2537                    
2538  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2539  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2540           if(nclouds_xz.ge.ncloxz_max)then           if(nclouds_xz.ge.ncloxz_max)then
2541              if(verbose)print*,              if(verbose.eq.1)print*,
2542       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2543       $           'XZ clouds exceeds vector dimention '       $           'XZ clouds exceeds vector dimention '
2544       $           ,'( ',ncloxz_max,' )'       $           ,'( ',ncloxz_max,' )'
2545  c     good2=.false.  c     good2=.false.
2546  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2547              do iv=1,nviews              do iv=1,nviews
2548                 mask_view(iv) = 6  c               mask_view(iv) = 6
2549                   mask_view(iv) =  mask_view(iv) + 2**5
2550              enddo              enddo
2551              iflag=1              iflag=1
2552              return              return
# Line 2208  c     goto 880         !fill ntp and go Line 2565  c     goto 880         !fill ntp and go
2565           enddo           enddo
2566           npt_tot=npt_tot+npt           npt_tot=npt_tot+npt
2567                    
2568           if(DEBUG)then           if(DEBUG.EQ.1)then
2569              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
2570              print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'                            print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'              
2571              print*,'- alfaxz1 ',alfaxz1_av(nclouds_xz)              print*,'- alfaxz1  ',alfaxz1_av(nclouds_xz)
2572              print*,'- alfaxz2 ',alfaxz2_av(nclouds_xz)              print*,'- alfaxz2  ',alfaxz2_av(nclouds_xz)
2573              print*,'- alfaxz3 ',alfaxz3_av(nclouds_xz)              print*,'- alfaxz3  ',alfaxz3_av(nclouds_xz)
2574              print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)              print*,'cp_useds1  ',(cp_useds1(icp),icp=1,ncp_tot)
2575              print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)              print*,'cp_useds2  ',(cp_useds2(icp),icp=1,ncp_tot)
2576                print*,'cpcloud_xz '
2577         $           ,(cpcloud_xz(nclouds_xz,icp),icp=1,ncp_tot)
2578              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)
2579  c$$$            print*,'nt-uple: ptcloud_xz(',nclouds_xz,') = '  c$$$            print*,'nt-uple: ptcloud_xz(',nclouds_xz,') = '
2580  c$$$     $           ,ptcloud_xz(nclouds_xz)  c$$$     $           ,ptcloud_xz(nclouds_xz)
# Line 2233  c$$$     $           ,(tr_cloud(iii),iii Line 2592  c$$$     $           ,(tr_cloud(iii),iii
2592           goto 91                           goto 91                
2593         endif                             endif                    
2594                
2595        if(DEBUG)then        if(DEBUG.EQ.1)then
2596           print*,'---------------------- '           print*,'---------------------- '
2597           print*,'X-Z total clouds ',nclouds_xz           print*,'X-Z total clouds ',nclouds_xz
2598           print*,' '           print*,' '
# Line 2254  c$$$     $           ,(tr_cloud(iii),iii Line 2613  c$$$     $           ,(tr_cloud(iii),iii
2613  **************************************************  **************************************************
2614    
2615        subroutine clouds_to_ctrack(iflag)        subroutine clouds_to_ctrack(iflag)
 c*****************************************************  
 c     02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
2616    
2617        include 'commontracker.f'        include 'commontracker.f'
2618        include 'level1.f'        include 'level1.f'
# Line 2264  c*************************************** Line 2620  c***************************************
2620        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2621        include 'common_mini_2.f'        include 'common_mini_2.f'
2622        include 'common_mech.f'        include 'common_mech.f'
2623  c      include 'momanhough_init.f'  
2624    
2625    
2626  *     output flag  *     output flag
# Line 2288  c      include 'momanhough_init.f' Line 2644  c      include 'momanhough_init.f'
2644  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2645  *     variables for track fitting  *     variables for track fitting
2646        double precision AL_INI(5)        double precision AL_INI(5)
 c      double precision tath  
2647  *     -----------------------------------------------------------  *     -----------------------------------------------------------
 c      real fitz(nplanes)        !z coordinates of the planes in cm  
2648    
2649          if(DEBUG.EQ.1)print*,'clouds_to_ctrack:'
2650    
2651    
2652        ntracks=0                 !counter of track candidates        ntracks=0                 !counter of track candidates
# Line 2313  c      real fitz(nplanes)        !z coor Line 2668  c      real fitz(nplanes)        !z coor
2668                 enddo                 enddo
2669              enddo              enddo
2670              ncp_ok=0              ncp_ok=0
2671              do icp=1,ncp_tot    !loop on couples              do icp=1,ncp_tot    !loop over couples
2672  *     get info on  *     get info on
2673                 cpintersec(icp)=min(                 cpintersec(icp)=min(
2674       $              cpcloud_yz(iyz,icp),       $              cpcloud_yz(iyz,icp),
# Line 2322  c      real fitz(nplanes)        !z coor Line 2677  c      real fitz(nplanes)        !z coor
2677       $    (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.
2678       $    (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.
2679       $              .false.)cpintersec(icp)=0       $              .false.)cpintersec(icp)=0
2680    *     cpintersec is >0 if yz and xz clouds contain the same image of couple icp
2681                 if(cpintersec(icp).ne.0)then                 if(cpintersec(icp).ne.0)then
2682                    ncp_ok=ncp_ok+1                      ncp_ok=ncp_ok+1  
2683                                        
# Line 2354  c      real fitz(nplanes)        !z coor Line 2710  c      real fitz(nplanes)        !z coor
2710                 nplused=nplused+ hit_plane(ip)                 nplused=nplused+ hit_plane(ip)
2711              enddo              enddo
2712                            
 c            if(nplused.lt.nplxz_min)goto 888 !next doublet  
             if(nplused.lt.nplyz_min)goto 888 !next doublet  
             if(ncp_ok.lt.ncpok)goto 888 !next cloud  
2713                            
2714              if(DEBUG)then              if(DEBUG.EQ.1)then
2715                 print*,'Combination ',iyz,ixz                 print*,'Combination ',iyz,ixz
2716       $              ,' db ',ptcloud_yz(iyz)       $              ,' db ',ptcloud_yz(iyz)
2717       $              ,' tr ',ptcloud_xz(ixz)       $              ,' tr ',ptcloud_xz(ixz)
2718       $              ,'  -----> # matching couples ',ncp_ok       $              ,'  -----> # matching couples ',ncp_ok
2719              endif              endif
2720    
2721    c            if(nplused.lt.nplxz_min)goto 888 !next combination
2722                if(nplused.lt.nplyz_min)goto 888 !next combination
2723                if(ncp_ok.lt.ncpok)goto 888 !next combination
2724    
2725  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
2726  c$$$  print*,'Configurazione cluster XZ'  c$$$  print*,'Configurazione cluster XZ'
2727  c$$$  print*,'1 -- ',(clx(1,i),i=1,ncp_plane(1))  c$$$  print*,'1 -- ',(clx(1,i),i=1,ncp_plane(1))
# Line 2409  c$$$            AL_INI(5) = (1.e2*alfaxz Line 2767  c$$$            AL_INI(5) = (1.e2*alfaxz
2767  c$$$              c$$$            
2768  c$$$            if(AL_INI(5).gt.defmax)goto 888 !next cloud  c$$$            if(AL_INI(5).gt.defmax)goto 888 !next cloud
2769                                                    
2770              if(DEBUG)then              if(DEBUG.EQ.1)then
2771                   print*,'track candidate', ntracks+1
2772                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))
2773                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))
2774                 print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4))                 print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4))
# Line 2442  c$$$            if(AL_INI(5).gt.defmax)g Line 2801  c$$$            if(AL_INI(5).gt.defmax)g
2801                                hit_plane(6)=icp6                                hit_plane(6)=icp6
2802                                if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6                                if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6
2803                                                                
2804                                  *                             ---------------------------------------
2805    *                             check if this group of couples has been
2806    *                             already fitted    
2807    *                             ---------------------------------------
2808                                  do ica=1,ntracks
2809                                     isthesame=1
2810                                     do ip=1,NPLANES
2811                                        if(hit_plane(ip).ne.0)then
2812                                           if(  CP_STORE(nplanes-ip+1,ica)
2813         $                                      .ne.
2814         $                                      cp_match(ip,hit_plane(ip)) )
2815         $                                      isthesame=0
2816                                        else
2817                                           if(  CP_STORE(nplanes-ip+1,ica)
2818         $                                      .ne.
2819         $                                      0 )
2820         $                                      isthesame=0
2821                                        endif
2822                                     enddo
2823                                     if(isthesame.eq.1)then
2824                                        if(DEBUG.eq.1)
2825         $                                   print*,'(already fitted)'
2826                                        goto 666 !jump to next combination
2827                                     endif
2828                                  enddo
2829    
2830                                call track_init !init TRACK common                                call track_init !init TRACK common
2831    
2832                                do ip=1,nplanes !loop on planes                                do ip=1,nplanes !loop on planes (bottom to top)
2833                                   if(hit_plane(ip).ne.0)then                                   if(hit_plane(ip).ne.0)then
2834                                      id=cp_match(ip,hit_plane(ip))                                      id=cp_match(ip,hit_plane(ip))
2835                                      is=is_cp(id)                                      is=is_cp(id)
# Line 2459  c$$$            if(AL_INI(5).gt.defmax)g Line 2843  c$$$            if(AL_INI(5).gt.defmax)g
2843  *                                   *************************  *                                   *************************
2844  c                                    call xyz_PAM(icx,icy,is,  c                                    call xyz_PAM(icx,icy,is,
2845  c     $                                   'COG2','COG2',0.,0.)  c     $                                   'COG2','COG2',0.,0.)
2846    c                                    call xyz_PAM(icx,icy,is, !(1)
2847    c     $                                   PFAdef,PFAdef,0.,0.) !(1)
2848                                      call xyz_PAM(icx,icy,is, !(1)                                      call xyz_PAM(icx,icy,is, !(1)
2849       $                                   PFAdef,PFAdef,0.,0.) !(1)       $                                   PFAdef,PFAdef,0.,0.,0.,0.)
2850  *                                   *************************  *                                   *************************
2851  *                                   -----------------------------  *                                   -----------------------------
2852                                      xgood(nplanes-ip+1)=1.                                      xgood(nplanes-ip+1)=1.
# Line 2487  c$$$                              enddo Line 2873  c$$$                              enddo
2873                                ifail=0 !error flag in chi^2 computation                                ifail=0 !error flag in chi^2 computation
2874                                jstep=0 !number of  minimization steps                                jstep=0 !number of  minimization steps
2875                                iprint=0                                iprint=0
2876  c                              if(DEBUG)iprint=1  c                              if(DEBUG.EQ.1)iprint=1
2877                                if(DEBUG)iprint=2                                if(DEBUG.EQ.1)iprint=2
2878                                call mini2(jstep,ifail,iprint)                                call mini2(jstep,ifail,iprint)
2879                                if(ifail.ne.0) then                                if(ifail.ne.0) then
2880                                   if(DEBUG)then                                   if(DEBUG.EQ.1)then
2881                                      print *,                                      print *,
2882       $                              '*** MINIMIZATION FAILURE *** '       $                              '*** MINIMIZATION FAILURE *** '
2883       $                              //'(clouds_to_ctrack)'       $                              //'(clouds_to_ctrack)'
# Line 2516  c                                 chi2=- Line 2902  c                                 chi2=-
2902  *     --------------------------  *     --------------------------
2903                                if(ntracks.eq.NTRACKSMAX)then                                if(ntracks.eq.NTRACKSMAX)then
2904                                                                    
2905                                   if(verbose)print*,                                   if(verbose.eq.1)print*,
2906       $                 '** warning ** number of candidate tracks '//       $                 '** warning ** number of candidate tracks '//
2907       $                 ' exceeds vector dimension '       $                 ' exceeds vector dimension '
2908       $                ,'( ',NTRACKSMAX,' )'       $                ,'( ',NTRACKSMAX,' )'
2909  c                                 good2=.false.  c                                 good2=.false.
2910  c                                 goto 880 !fill ntp and go to next event                      c                                 goto 880 !fill ntp and go to next event                    
2911                                   do iv=1,nviews                                   do iv=1,nviews
2912                                      mask_view(iv) = 7  c                                    mask_view(iv) = 7
2913                                        mask_view(iv) = mask_view(iv) + 2**6
2914                                   enddo                                   enddo
2915                                   iflag=1                                   iflag=1
2916                                   return                                   return
# Line 2531  c                                 goto 8 Line 2918  c                                 goto 8
2918                                                                
2919                                ntracks = ntracks + 1                                ntracks = ntracks + 1
2920                                                                
2921  c$$$                              ndof=0                                                                do ip=1,nplanes !top to bottom
2922                                do ip=1,nplanes  
 c$$$                                 ndof=ndof  
 c$$$     $                                +int(xgood(ip))  
 c$$$     $                                +int(ygood(ip))  
2923                                   XV_STORE(ip,ntracks)=sngl(xv(ip))                                   XV_STORE(ip,ntracks)=sngl(xv(ip))
2924                                   YV_STORE(ip,ntracks)=sngl(yv(ip))                                   YV_STORE(ip,ntracks)=sngl(yv(ip))
2925                                   ZV_STORE(ip,ntracks)=sngl(zv(ip))                                                                       ZV_STORE(ip,ntracks)=sngl(zv(ip))                                    
# Line 2551  c$$$     $                               Line 2935  c$$$     $                              
2935                                   AYV_STORE(ip,ntracks)=sngl(ayv(ip))                                   AYV_STORE(ip,ntracks)=sngl(ayv(ip))
2936                                   XGOOD_STORE(ip,ntracks)=sngl(xgood(ip))                                   XGOOD_STORE(ip,ntracks)=sngl(xgood(ip))
2937                                   YGOOD_STORE(ip,ntracks)=sngl(ygood(ip))                                   YGOOD_STORE(ip,ntracks)=sngl(ygood(ip))
2938    *                                NB! hit_plane is defined from bottom to top
2939                                   if(hit_plane(ip).ne.0)then                                   if(hit_plane(ip).ne.0)then
2940                                      CP_STORE(nplanes-ip+1,ntracks)=                                      CP_STORE(nplanes-ip+1,ntracks)=
2941       $                                   cp_match(ip,hit_plane(ip))       $                                   cp_match(ip,hit_plane(ip))
2942                                        SENSOR_STORE(nplanes-ip+1,ntracks)
2943         $                              = is_cp(cp_match(ip,hit_plane(ip)))
2944                                        LADDER_STORE(nplanes-ip+1,ntracks)
2945         $                                   = LADDER(
2946         $                                   clx(ip,icp_cp(
2947         $                                   cp_match(ip,hit_plane(ip)
2948         $                                   ))));
2949                                   else                                   else
2950                                      CP_STORE(nplanes-ip+1,ntracks)=0                                      CP_STORE(nplanes-ip+1,ntracks)=0
2951                                        SENSOR_STORE(nplanes-ip+1,ntracks)=0
2952                                        LADDER_STORE(nplanes-ip+1,ntracks)=0
2953                                   endif                                   endif
2954                                   CLS_STORE(nplanes-ip+1,ntracks)=0                                   BX_STORE(ip,ntracks)=0!I dont need it now
2955                                     BY_STORE(ip,ntracks)=0!I dont need it now
2956                                     CLS_STORE(ip,ntracks)=0
2957                                   do i=1,5                                   do i=1,5
2958                                      AL_STORE(i,ntracks)=sngl(AL(i))                                      AL_STORE(i,ntracks)=sngl(AL(i))
2959                                   enddo                                   enddo
2960                                enddo                                enddo
2961                                                                
 c$$$  *                             Number of Degree Of Freedom  
 c$$$  ndof=ndof-5                            
 c$$$  *                             reduced chi^2  
 c$$$  rchi2=chi2/dble(ndof)  
2962                                RCHI2_STORE(ntracks)=chi2                                RCHI2_STORE(ntracks)=chi2
2963                                                                
2964  *     --------------------------------  *     --------------------------------
# Line 2590  c$$$  rchi2=chi2/dble(ndof) Line 2982  c$$$  rchi2=chi2/dble(ndof)
2982           return           return
2983        endif        endif
2984                
2985        if(DEBUG)then  c$$$      if(DEBUG.EQ.1)then
2986           print*,'****** TRACK CANDIDATES ***********'  c$$$         print*,'****** TRACK CANDIDATES ***********'
2987           print*,'#         R. chi2        RIG'  c$$$         print*,'#         R. chi2        RIG'
2988           do i=1,ntracks  c$$$         do i=1,ntracks
2989              print*,i,' --- ',rchi2_store(i),' --- '  c$$$            print*,i,' --- ',rchi2_store(i),' --- '
2990       $           ,1./abs(AL_STORE(5,i))  c$$$     $           ,1./abs(AL_STORE(5,i))
2991           enddo  c$$$         enddo
2992           print*,'***********************************'  c$$$         print*,'***********************************'
2993    c$$$      endif
2994          if(DEBUG.EQ.1)then
2995            print*,'****** TRACK CANDIDATES *****************'
2996            print*,'#         R. chi2        RIG         ndof'
2997            do i=1,ntracks
2998              ndof=0                !(1)
2999              do ii=1,nplanes       !(1)
3000                ndof=ndof           !(1)
3001         $           +int(xgood_store(ii,i)) !(1)
3002         $           +int(ygood_store(ii,i)) !(1)
3003              enddo                 !(1)
3004              print*,i,' --- ',rchi2_store(i),' --- '
3005         $         ,1./abs(AL_STORE(5,i)),' --- ',ndof
3006            enddo
3007            print*,'*****************************************'
3008        endif        endif
3009                
3010                
# Line 2616  c$$$  rchi2=chi2/dble(ndof) Line 3023  c$$$  rchi2=chi2/dble(ndof)
3023    
3024        subroutine refine_track(ibest)        subroutine refine_track(ibest)
3025    
 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******************************************************  
3026    
3027        include 'commontracker.f'        include 'commontracker.f'
3028        include 'level1.f'        include 'level1.f'
# Line 2628  c*************************************** Line 3030  c***************************************
3030        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3031        include 'common_mini_2.f'        include 'common_mini_2.f'
3032        include 'common_mech.f'        include 'common_mech.f'
 c      include 'momanhough_init.f'  
 c      include 'level1.f'  
3033        include 'calib.f'        include 'calib.f'
3034    
   
3035  *     flag to chose PFA  *     flag to chose PFA
3036        character*10 PFA        character*10 PFA
3037        common/FINALPFA/PFA        common/FINALPFA/PFA
3038    
3039          real k(6)
3040          DATA k/1.099730,0.418900,0.220939,0.220907,0.418771,1.100674/
3041    
3042          real xp,yp,zp
3043          real xyzp(3),bxyz(3)
3044          equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3))
3045    
3046          if(DEBUG.EQ.1)print*,'refine_track:'
3047  *     =================================================  *     =================================================
3048  *     new estimate of positions using ETA algorithm  *     new estimate of positions using ETA algorithm
3049  *                          and  *                          and
# Line 2645  c      include 'level1.f' Line 3052  c      include 'level1.f'
3052        call track_init        call track_init
3053        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3054    
3055             xP=XV_STORE(nplanes-ip+1,ibest)
3056             yP=YV_STORE(nplanes-ip+1,ibest)
3057             zP=ZV_STORE(nplanes-ip+1,ibest)
3058             call gufld(xyzp,bxyz)
3059             BX_STORE(nplanes-ip+1,ibest)=bxyz(1)
3060             BY_STORE(nplanes-ip+1,ibest)=bxyz(2)
3061    c$$$  bxyz(1)=0
3062    c$$$         bxyz(2)=0
3063    c$$$         bxyz(3)=0
3064  *     |||||||||||||||||||||||||||||||||||||||||||||||||  *     |||||||||||||||||||||||||||||||||||||||||||||||||
3065  *     -------------------------------------------------  *     -------------------------------------------------
3066  *     If the plane has been already included, it just  *     If the plane has been already included, it just
# Line 2665  c      include 'level1.f' Line 3081  c      include 'level1.f'
3081       $           ,ip_cp(id),ip       $           ,ip_cp(id),ip
3082              icx=clx(ip,icp)              icx=clx(ip,icp)
3083              icy=cly(ip,icp)              icy=cly(ip,icp)
3084    c            call xyz_PAM(icx,icy,is,
3085    c     $           PFA,PFA,
3086    c     $           AXV_STORE(nplanes-ip+1,ibest),
3087    c     $           AYV_STORE(nplanes-ip+1,ibest))
3088              call xyz_PAM(icx,icy,is,              call xyz_PAM(icx,icy,is,
 c     $           'ETA2','ETA2',  
3089       $           PFA,PFA,       $           PFA,PFA,
3090       $           AXV_STORE(nplanes-ip+1,ibest),       $           AXV_STORE(nplanes-ip+1,ibest),
3091       $           AYV_STORE(nplanes-ip+1,ibest))       $           AYV_STORE(nplanes-ip+1,ibest),
3092  c$$$  call xyz_PAM(icx,icy,is,       $           bxyz(1),
3093  c$$$  $              'COG2','COG2',       $           bxyz(2)
3094  c$$$  $              0.,       $           )
3095  c$$$  $              0.)  
3096              xm(nplanes-ip+1) = xPAM              xm(nplanes-ip+1) = xPAM
3097              ym(nplanes-ip+1) = yPAM              ym(nplanes-ip+1) = yPAM
3098              zm(nplanes-ip+1) = zPAM              zm(nplanes-ip+1) = zPAM
# Line 2682  c$$$  $              0.) Line 3101  c$$$  $              0.)
3101              resx(nplanes-ip+1) = resxPAM              resx(nplanes-ip+1) = resxPAM
3102              resy(nplanes-ip+1) = resyPAM              resy(nplanes-ip+1) = resyPAM
3103    
3104  c            dedxtrk(nplanes-ip+1) = (dedx(icx)+dedx(icy))/2. !(1)              dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx))
3105              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)  
3106                            
3107  *     |||||||||||||||||||||||||||||||||||||||||||||||||  *     |||||||||||||||||||||||||||||||||||||||||||||||||
3108  *     -------------------------------------------------  *     -------------------------------------------------
# Line 2699  c            dedxtrk(nplanes-ip+1) = (de Line 3117  c            dedxtrk(nplanes-ip+1) = (de
3117                                
3118  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3119  *     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)  
3120              call whichsensor(ip,xP,yP,nldt,ist)              call whichsensor(ip,xP,yP,nldt,ist)
3121  *     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
3122              if(nldt.eq.0.or.ist.eq.0)goto 133              if(nldt.eq.0.or.ist.eq.0)goto 133
3123    
3124                SENSOR_STORE(nplanes-ip+1,IBEST)=ist
3125                LADDER_STORE(nplanes-ip+1,IBEST)=nldt
3126  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3127    
3128              if(DEBUG)then              if(DEBUG.EQ.1)then
3129                 print*,                 print*,
3130       $              '------ Plane ',ip,' intersected on LADDER ',nldt       $              '------ Plane ',ip,' intersected on LADDER ',nldt
3131       $              ,' SENSOR ',ist       $              ,' SENSOR ',ist
# Line 2718  c            dedxtrk(nplanes-ip+1) = (de Line 3136  c            dedxtrk(nplanes-ip+1) = (de
3136  *     ===========================================  *     ===========================================
3137  *     STEP 1 >>>>>>>  try to include a new couple  *     STEP 1 >>>>>>>  try to include a new couple
3138  *     ===========================================  *     ===========================================
3139  c            if(DEBUG)print*,'>>>> try to include a new couple'  c            if(DEBUG.EQ.1)print*,'>>>> try to include a new couple'
3140              distmin=1000000.              distmin=1000000.
3141              xmm = 0.              xmm = 0.
3142              ymm = 0.              ymm = 0.
# Line 2741  c     $              cl_used(icy).eq.1.o Line 3159  c     $              cl_used(icy).eq.1.o
3159  *            *          
3160                 call xyz_PAM(icx,icy,ist,                 call xyz_PAM(icx,icy,ist,
3161       $              PFA,PFA,       $              PFA,PFA,
 c     $              'ETA2','ETA2',  
3162       $              AXV_STORE(nplanes-ip+1,ibest),       $              AXV_STORE(nplanes-ip+1,ibest),
3163       $              AYV_STORE(nplanes-ip+1,ibest))       $              AYV_STORE(nplanes-ip+1,ibest),
3164         $              bxyz(1),
3165         $              bxyz(2)
3166         $              )
3167                                
3168                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3169                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3170                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
3171                 if(DEBUG)print*,'( couple ',id                 if(DEBUG.EQ.1)print*,'( couple ',id
3172       $              ,' ) normalized distance ',distance       $              ,' ) distance ',distance
3173                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3174                    xmm = xPAM                    xmm = xPAM
3175                    ymm = yPAM                    ymm = yPAM
# Line 2758  c     $              'ETA2','ETA2', Line 3178  c     $              'ETA2','ETA2',
3178                    rymm = resyPAM                    rymm = resyPAM
3179                    distmin = distance                    distmin = distance
3180                    idm = id                                      idm = id                  
3181  c                 dedxmm = (dedx(icx)+dedx(icy))/2. !(1)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3182                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3183                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)  c     QUIQUI --> non devo moltiplicare per clinc?!?!?!
3184                      clincnewc=10*sqrt(rymm**2+rxmm**2 !QUIQUI
3185         $                 +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2))) !QUIQUI
3186                 endif                 endif
3187   1188          continue   1188          continue
3188              enddo               !end loop on couples on plane icp              enddo               !end loop on couples on plane icp
3189              if(distmin.le.clinc)then                    c            if(distmin.le.clinc)then     !QUIQUI              
3190                if(distmin.le.clincnewc)then     !QUIQUI              
3191  *              -----------------------------------  *              -----------------------------------
3192                 xm(nplanes-ip+1) = xmm         !<<<                 xm(nplanes-ip+1) = xmm !<<<
3193                 ym(nplanes-ip+1) = ymm         !<<<                 ym(nplanes-ip+1) = ymm !<<<
3194                 zm(nplanes-ip+1) = zmm         !<<<                 zm(nplanes-ip+1) = zmm !<<<
3195                 xgood(nplanes-ip+1) = 1        !<<<                 xgood(nplanes-ip+1) = 1 !<<<
3196                 ygood(nplanes-ip+1) = 1        !<<<                 ygood(nplanes-ip+1) = 1 !<<<
3197                 resx(nplanes-ip+1)=rxmm        !<<<                 resx(nplanes-ip+1)=rxmm !<<<
3198                 resy(nplanes-ip+1)=rymm        !<<<                 resy(nplanes-ip+1)=rymm !<<<
3199  c              dedxtrk(nplanes-ip+1) = dedxmm !<<<  !(1)                 dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3200                 dedxtrk_x(nplanes-ip+1) = dedxmmx    !(1)                 dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
                dedxtrk_y(nplanes-ip+1) = dedxmmy    !(1)  
3201  *              -----------------------------------  *              -----------------------------------
3202                 CP_STORE(nplanes-ip+1,ibest)=idm                       CP_STORE(nplanes-ip+1,ibest)=idm      
3203                 if(DEBUG)print*,'%%%% included couple ',idm                 if(DEBUG.EQ.1)print*,'%%%% included couple ',idm
3204       $              ,' (norm.dist.= ',distmin,', cut ',clinc,' )'       $              ,' (dist.= ',distmin,', cut ',clinc,' )'
3205                 goto 133         !next plane                 goto 133         !next plane
3206              endif              endif
3207  *     ================================================  *     ================================================
3208  *     STEP 2 >>>>>>>  try to include a single cluster  *     STEP 2 >>>>>>>  try to include a single cluster
3209  *                     either from a couple or single  *                     either from a couple or single
3210  *     ================================================  *     ================================================
3211  c            if(DEBUG)print*,'>>>> try to include a new cluster'  c            if(DEBUG.EQ.1)print*,'>>>> try to include a new cluster'
3212              distmin=1000000.              distmin=1000000.
3213              xmm_A = 0.          !---------------------------              xmm_A = 0.          !---------------------------
3214              ymm_A = 0.          ! init variables that              ymm_A = 0.          ! init variables that
# Line 2812  c            if(DEBUG)print*,'>>>> try t Line 3234  c            if(DEBUG)print*,'>>>> try t
3234  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
3235                 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)
3236  *                                              !jump to the Y cluster  *                                              !jump to the Y cluster
3237    c               call xyz_PAM(icx,0,ist,
3238    c     $              PFA,PFA,
3239    c     $              AXV_STORE(nplanes-ip+1,ibest),0.)              
3240                 call xyz_PAM(icx,0,ist,                 call xyz_PAM(icx,0,ist,
 c     $              'ETA2','ETA2',  
3241       $              PFA,PFA,       $              PFA,PFA,
3242       $              AXV_STORE(nplanes-ip+1,ibest),0.)                     $              AXV_STORE(nplanes-ip+1,ibest),0.,
3243         $              bxyz(1),
3244         $              bxyz(2)
3245         $              )              
3246                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3247                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3248                 if(DEBUG)print*,'( cl-X ',icx                 if(DEBUG.EQ.1)print*,'( cl-X ',icx
3249       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3250                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3251                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3252                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 2831  c     $              'ETA2','ETA2', Line 3258  c     $              'ETA2','ETA2',
3258                    rymm = resyPAM                    rymm = resyPAM
3259                    distmin = distance                    distmin = distance
3260                    iclm = icx                    iclm = icx
3261  c                  dedxmm = dedx(icx) !(1)  c                  dedxmm = sgnl(icx) !(1)
3262                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3263                    dedxmmy = 0.        !(1)                    dedxmmy = 0.        !(1)
3264                 endif                                   endif                  
3265  11881          continue  11881          continue
# Line 2840  c                  dedxmm = dedx(icx) !( Line 3267  c                  dedxmm = dedx(icx) !(
3267  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
3268                 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)
3269  *                                              !jump to the next couple  *                                              !jump to the next couple
3270    c               call xyz_PAM(0,icy,ist,
3271    c     $              PFA,PFA,
3272    c     $              0.,AYV_STORE(nplanes-ip+1,ibest))
3273                 call xyz_PAM(0,icy,ist,                 call xyz_PAM(0,icy,ist,
 c     $              'ETA2','ETA2',  
3274       $              PFA,PFA,       $              PFA,PFA,
3275       $              0.,AYV_STORE(nplanes-ip+1,ibest))       $              0.,AYV_STORE(nplanes-ip+1,ibest),
3276         $              bxyz(1),
3277         $              bxyz(2)
3278         $              )
3279                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3280                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3281                 if(DEBUG)print*,'( cl-Y ',icy                 if(DEBUG.EQ.1)print*,'( cl-Y ',icy
3282       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3283                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3284                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3285                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 2859  c     $              'ETA2','ETA2', Line 3291  c     $              'ETA2','ETA2',
3291                    rymm = resyPAM                    rymm = resyPAM
3292                    distmin = distance                    distmin = distance
3293                    iclm = icy                    iclm = icy
3294  c                 dedxmm = dedx(icy)  !(1)  c                 dedxmm = sgnl(icy)  !(1)
3295                    dedxmmx = 0.        !(1)                    dedxmmx = 0.        !(1)
3296                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3297                 endif                                   endif                  
3298  11882          continue  11882          continue
3299              enddo               !end loop on cluster inside couples              enddo               !end loop on cluster inside couples
3300  *----- single clusters -----------------------------------------------      *----- single clusters -----------------------------------------------  
3301    c            print*,'## ncls(',ip,') ',ncls(ip)
3302              do ic=1,ncls(ip)    !loop on single clusters              do ic=1,ncls(ip)    !loop on single clusters
3303                 icl=cls(ip,ic)                 icl=cls(ip,ic)
3304  c               if(cl_used(icl).eq.1.or.     !if the cluster is already used  c               if(cl_used(icl).eq.1.or.     !if the cluster is already used
# Line 2874  c               if(cl_used(icl).eq.1.or. Line 3307  c               if(cl_used(icl).eq.1.or.
3307       $              .false.)goto 18882      !jump to the next singlet       $              .false.)goto 18882      !jump to the next singlet
3308                 if(mod(VIEW(icl),2).eq.0)then!<---- X view                 if(mod(VIEW(icl),2).eq.0)then!<---- X view
3309                    call xyz_PAM(icl,0,ist,                    call xyz_PAM(icl,0,ist,
 c     $                 'ETA2','ETA2',  
3310       $                 PFA,PFA,       $                 PFA,PFA,
3311       $                 AXV_STORE(nplanes-ip+1,ibest),0.)       $                 AXV_STORE(nplanes-ip+1,ibest),0.,
3312         $                 bxyz(1),
3313         $                 bxyz(2)
3314         $                 )
3315                 else                         !<---- Y view                 else                         !<---- Y view
3316                    call xyz_PAM(0,icl,ist,                    call xyz_PAM(0,icl,ist,
 c     $                 'ETA2','ETA2',  
3317       $                 PFA,PFA,       $                 PFA,PFA,
3318       $                 0.,AYV_STORE(nplanes-ip+1,ibest))       $                 0.,AYV_STORE(nplanes-ip+1,ibest),
3319         $                 bxyz(1),
3320         $                 bxyz(2)
3321         $                 )
3322                 endif                 endif
3323    
3324                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3325                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3326                 if(DEBUG)print*,'( cl-s ',icl                 if(DEBUG.EQ.1)print*,'( cl-s ',icl
3327       $              ,' ) normalized distance ',distance       $              ,' ) distance ',distance
3328                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3329    c                  if(DEBUG.EQ.1)print*,'YES'
3330                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3331                    ymm_A = yPAM_A                    ymm_A = yPAM_A
3332                    zmm_A = zPAM_A                    zmm_A = zPAM_A
# Line 2899  c     $                 'ETA2','ETA2', Line 3337  c     $                 'ETA2','ETA2',
3337                    rymm = resyPAM                    rymm = resyPAM
3338                    distmin = distance                      distmin = distance  
3339                    iclm = icl                    iclm = icl
 c                  dedxmm = dedx(icl)                   !(1)  
3340                    if(mod(VIEW(icl),2).eq.0)then !<---- X view                    if(mod(VIEW(icl),2).eq.0)then !<---- X view
3341                       dedxmmx = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmx = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3342                       dedxmmy = 0.                       !(1)                       dedxmmy = 0.                  
3343                    else          !<---- Y view                    else          !<---- Y view
3344                       dedxmmx = 0.                       !(1)                       dedxmmx = 0.                  
3345                       dedxmmy = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmy = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3346                    endif                    endif
3347                 endif                                   endif                  
3348  18882          continue  18882          continue
3349              enddo               !end loop on single clusters              enddo               !end loop on single clusters
3350    c            print*,'## distmin ', distmin,' clinc ',clinc
3351    
3352              if(distmin.le.clinc)then                    c     QUIQUI------------
3353                  c     anche qui: non ci vuole clinc???
3354                 CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<                  if(iclm.ne.0)then
 *              ----------------------------  
 c               print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'  
3355                 if(mod(VIEW(iclm),2).eq.0)then                 if(mod(VIEW(iclm),2).eq.0)then
3356                    XGOOD(nplanes-ip+1)=1.                    clincnew=
3357                    resx(nplanes-ip+1)=rxmm       $                 20*
3358                    if(DEBUG)print*,'%%%% included X-cl ',iclm       $                 sqrt(rxmm**2+RCHI2_STORE(ibest)*k(ip)*cov(1,1))
3359  c                  if(.true.)print*,'%%%% included X-cl ',iclm                 else if(mod(VIEW(iclm),2).ne.0)then
3360       $                 ,'( chi^2, ',RCHI2_STORE(ibest)                    clincnew=
3361       $                 ,', norm.dist.= ',distmin       $                 10*
3362       $                 ,', cut ',clinc,' )'       $                 sqrt(rymm**2+RCHI2_STORE(ibest)*k(ip)*cov(2,2))
3363                 else                 endif
3364                    YGOOD(nplanes-ip+1)=1.  c     QUIQUI------------
3365                    resy(nplanes-ip+1)=rymm                
3366                    if(DEBUG)print*,'%%%% included Y-cl ',iclm                 if(distmin.le.clincnew)then   !QUIQUI
3367  c                  if(.true.)print*,'%%%% included Y-cl ',iclm  c     if(distmin.le.clinc)then          !QUIQUI          
3368       $                 ,'( chi^2, ',RCHI2_STORE(ibest)                    
3369       $                 ,', norm.dist.= ', distmin                    CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<    
3370       $                 ,', cut ',clinc,' )'  *     ----------------------------
3371    c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
3372                      if(mod(VIEW(iclm),2).eq.0)then
3373                         XGOOD(nplanes-ip+1)=1.
3374                         resx(nplanes-ip+1)=rxmm
3375                         if(DEBUG.EQ.1)print*,'%%%% included X-cl ',iclm
3376         $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3377         $                    ,', dist.= ',distmin
3378         $                    ,', cut ',clinc,' )'
3379                      else
3380                         YGOOD(nplanes-ip+1)=1.
3381                         resy(nplanes-ip+1)=rymm
3382                         if(DEBUG.EQ.1)print*,'%%%% included Y-cl ',iclm
3383         $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3384         $                    ,', dist.= ', distmin
3385         $                    ,', cut ',clinc,' )'
3386                      endif
3387    c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
3388    *     ----------------------------
3389                      xm_A(nplanes-ip+1) = xmm_A
3390                      ym_A(nplanes-ip+1) = ymm_A
3391                      xm_B(nplanes-ip+1) = xmm_B
3392                      ym_B(nplanes-ip+1) = ymm_B
3393                      zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.
3394                      dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3395                      dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
3396    *     ----------------------------
3397                 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)  
 *              ----------------------------  
3398              endif              endif
3399           endif           endif
3400   133     continue   133     continue
# Line 2962  c              dedxtrk(nplanes-ip+1) = d Line 3413  c              dedxtrk(nplanes-ip+1) = d
3413  *                                                 *  *                                                 *
3414  *                                                 *  *                                                 *
3415  **************************************************  **************************************************
 cccccc 12/08/2006 modified by elena ---> (1)  
3416  *  *
3417        subroutine clean_XYclouds(ibest,iflag)        subroutine clean_XYclouds(ibest,iflag)
3418    
3419        include 'commontracker.f'        include 'commontracker.f'
3420        include 'level1.f'        include 'level1.f'
3421        include 'common_momanhough.f'        include 'common_momanhough.f'
3422  c      include 'momanhough_init.f'        include 'level2.f'      
       include 'level2.f'        !(1)  
 c      include 'calib.f'  
 c      include 'level1.f'  
   
3423    
3424          if(DEBUG.EQ.1)print*,'clean_XYclouds:'
3425    
3426        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3427    
# Line 2984  c      include 'level1.f' Line 3431  c      include 'level1.f'
3431              if(id.ne.0)then              if(id.ne.0)then
3432                 iclx=clx(ip,icp_cp(id))                 iclx=clx(ip,icp_cp(id))
3433                 icly=cly(ip,icp_cp(id))                 icly=cly(ip,icp_cp(id))
3434  c               cl_used(iclx)=1  !tag used clusters  c$$$               cl_used(iclx)=ntrk  !tag used clusters
3435  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)  
3436              elseif(icl.ne.0)then              elseif(icl.ne.0)then
3437  c               cl_used(icl)=1   !tag used clusters  c$$$               cl_used(icl)=ntrk   !tag used clusters
                cl_used(icl)=ntrk   !tag used clusters !1)  
3438              endif              endif
3439                            
 c               if(DEBUG)then  
 c                  print*,ip,' <<< ',id  
 c               endif  
3440  *     -----------------------------  *     -----------------------------
3441  *     remove the couple from clouds  *     remove the couple from clouds
3442  *     remove also vitual couples containing the  *     remove also vitual couples containing the
# Line 3012  c               endif Line 3453  c               endif
3453       $              cly(ip,icp).eq.icl       $              cly(ip,icp).eq.icl
3454       $              )then       $              )then
3455                    id=id_cp(ip,icp,1)                    id=id_cp(ip,icp,1)
3456                    if(DEBUG)then                    if(DEBUG.EQ.1)then
3457                       print*,ip,' <<< cp ',id                       print*,ip,' <<< cp ',id
3458       $                    ,' ( cl-x '       $                    ,' ( cl-x '
3459       $                    ,clx(ip,icp)       $                    ,clx(ip,icp)
# Line 3056  c               endif Line 3497  c               endif
3497        include 'level1.f'        include 'level1.f'
3498        include 'common_momanhough.f'        include 'common_momanhough.f'
3499        include 'level2.f'        include 'level2.f'
 c      include 'level1.f'  
3500    
3501    *     ---------------------------------
3502    *     variables initialized from level1
3503    *     ---------------------------------
3504        do i=1,nviews        do i=1,nviews
3505           good2(i)=good1(i)           good2(i)=good1(i)
3506             do j=1,nva1_view
3507                vkflag(i,j)=1
3508                if(cnnev(i,j).le.0)then
3509                   vkflag(i,j)=cnnev(i,j)
3510                endif
3511             enddo
3512        enddo        enddo
3513    *     ----------------
3514    *     level2 variables
3515    *     ----------------
3516        NTRK = 0        NTRK = 0
3517        do it=1,NTRKMAX        do it=1,NTRKMAX
3518           IMAGE(IT)=0           IMAGE(IT)=0
# Line 3073  c      include 'level1.f' Line 3523  c      include 'level1.f'
3523              ZM_nt(IP,IT) = 0              ZM_nt(IP,IT) = 0
3524              RESX_nt(IP,IT) = 0              RESX_nt(IP,IT) = 0
3525              RESY_nt(IP,IT) = 0              RESY_nt(IP,IT) = 0
3526                TAILX_nt(IP,IT) = 0
3527                TAILY_nt(IP,IT) = 0
3528                XBAD(IP,IT) = 0
3529                YBAD(IP,IT) = 0
3530              XGOOD_nt(IP,IT) = 0              XGOOD_nt(IP,IT) = 0
3531              YGOOD_nt(IP,IT) = 0              YGOOD_nt(IP,IT) = 0
3532                LS(IP,IT) = 0
3533              DEDX_X(IP,IT) = 0              DEDX_X(IP,IT) = 0
3534              DEDX_Y(IP,IT) = 0              DEDX_Y(IP,IT) = 0
3535              CLTRX(IP,IT) = 0              CLTRX(IP,IT) = 0
3536              CLTRY(IP,IT) = 0              CLTRY(IP,IT) = 0
3537                multmaxx(ip,it) = 0
3538                seedx(ip,it)    = 0
3539                xpu(ip,it)      = 0
3540                multmaxy(ip,it) = 0
3541                seedy(ip,it)    = 0
3542                ypu(ip,it)      = 0
3543           enddo           enddo
3544           do ipa=1,5           do ipa=1,5
3545              AL_nt(IPA,IT) = 0              AL_nt(IPA,IT) = 0
# Line 3207  c      include 'level1.f' Line 3668  c      include 'level1.f'
3668    
3669            
3670        include 'commontracker.f'        include 'commontracker.f'
 c      include 'level1.f'  
3671        include 'level1.f'        include 'level1.f'
3672        include 'common_momanhough.f'        include 'common_momanhough.f'
3673        include 'level2.f'        include 'level2.f'
3674        include 'common_mini_2.f'        include 'common_mini_2.f'
3675        real sinth,phi,pig              include 'calib.f'
3676    
3677          character*10 PFA
3678          common/FINALPFA/PFA
3679    
3680          real sinth,phi,pig
3681          integer ssensor,sladder
3682        pig=acos(-1.)        pig=acos(-1.)
3683    
3684    *     -------------------------------------
3685        chi2_nt(ntr)        = sngl(chi2)        chi2_nt(ntr)        = sngl(chi2)
3686        nstep_nt(ntr)       = nstep        nstep_nt(ntr)       = nstep
3687    *     -------------------------------------
3688        phi   = al(4)                  phi   = al(4)          
3689        sinth = al(3)                    sinth = al(3)            
3690        if(sinth.lt.0)then              if(sinth.lt.0)then      
# Line 3230  c      include 'level1.f' Line 3697  c      include 'level1.f'
3697       $     phi = phi + 2*pig         $     phi = phi + 2*pig  
3698        al(4) = phi                      al(4) = phi              
3699        al(3) = sinth                    al(3) = sinth            
   
3700        do i=1,5        do i=1,5
3701           al_nt(i,ntr)     = sngl(al(i))           al_nt(i,ntr)     = sngl(al(i))
3702           do j=1,5           do j=1,5
3703              coval(i,j,ntr) = sngl(cov(i,j))              coval(i,j,ntr) = sngl(cov(i,j))
3704           enddo           enddo
3705        enddo        enddo
3706          *     -------------------------------------      
3707        do ip=1,nplanes           ! loop on planes        do ip=1,nplanes           ! loop on planes
3708           xgood_nt(ip,ntr) = int(xgood(ip))           xgood_nt(ip,ntr) = int(xgood(ip))
3709           ygood_nt(ip,ntr) = int(ygood(ip))           ygood_nt(ip,ntr) = int(ygood(ip))
# Line 3246  c      include 'level1.f' Line 3712  c      include 'level1.f'
3712           zm_nt(ip,ntr)    = sngl(zm(ip))           zm_nt(ip,ntr)    = sngl(zm(ip))
3713           RESX_nt(IP,ntr)  = sngl(resx(ip))           RESX_nt(IP,ntr)  = sngl(resx(ip))
3714           RESY_nt(IP,ntr)  = sngl(resy(ip))           RESY_nt(IP,ntr)  = sngl(resy(ip))
3715             TAILX_nt(IP,ntr) = 0.
3716             TAILY_nt(IP,ntr) = 0.
3717           xv_nt(ip,ntr)    = sngl(xv(ip))           xv_nt(ip,ntr)    = sngl(xv(ip))
3718           yv_nt(ip,ntr)    = sngl(yv(ip))           yv_nt(ip,ntr)    = sngl(yv(ip))
3719           zv_nt(ip,ntr)    = sngl(zv(ip))           zv_nt(ip,ntr)    = sngl(zv(ip))
3720           axv_nt(ip,ntr)   = sngl(axv(ip))           axv_nt(ip,ntr)   = sngl(axv(ip))
3721           ayv_nt(ip,ntr)   = sngl(ayv(ip))           ayv_nt(ip,ntr)   = sngl(ayv(ip))  
3722           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)  
3723           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)             factor = sqrt(
3724         $        tan( acos(-1.) * sngl(axv(ip)) /180. )**2 +
3725         $        tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 +
3726         $        1. )
3727    
3728             dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)
3729             dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)  
3730        
3731           id  = CP_STORE(ip,IDCAND)           ax   = axv_nt(ip,ntr)
3732             ay   = ayv_nt(ip,ntr)
3733             bfx  = BX_STORE(ip,IDCAND)
3734             bfy  = BY_STORE(ip,IDCAND)
3735    c$$$         if(ip.eq.6) ax = -1. * axv_nt(ip,ntr)
3736    c$$$         if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND)
3737    c$$$         tgtemp   = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001
3738    c$$$         angx     = 180.*atan(tgtemp)/acos(-1.)
3739    c$$$         tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
3740    c$$$         angy    = 180.*atan(tgtemp)/acos(-1.)
3741    
3742             angx = effectiveangle(ax,2*ip,bfy)
3743             angy = effectiveangle(ay,2*ip-1,bfx)
3744            
3745            
3746    c         print*,'* ',ip,bfx,bfy,angx,angy
3747    
3748             id  = CP_STORE(ip,IDCAND) ! couple id
3749           icl = CLS_STORE(ip,IDCAND)           icl = CLS_STORE(ip,IDCAND)
3750             ssensor = -1
3751             sladder = -1
3752             ssensor = SENSOR_STORE(ip,IDCAND)
3753             sladder = LADDER_STORE(ip,IDCAND)
3754             if(ip.eq.6.and.ssensor.ne.0)ssensor = 3 - ssensor !notazione paolo x align
3755             LS(IP,ntr)      = ssensor+10*sladder
3756    
3757           if(id.ne.0)then           if(id.ne.0)then
3758    c           >>> is a couple
3759              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))
3760              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))
3761  c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)  
3762                cl_used(cltrx(ip,ntr)) = 1     !tag used clusters          
3763                cl_used(cltry(ip,ntr)) = 1     !tag used clusters          
3764    
3765                xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))
3766                ybad(ip,ntr)= nbadstrips(4,cly(nplanes-ip+1,icp_cp(id)))
3767    
3768    
3769                if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0)
3770         $           dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3771                if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)
3772         $           dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3773    
3774                multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
3775         $                         +10000*mult(cltrx(ip,ntr))
3776                seedx(ip,ntr)    = clsignal(indmax(cltrx(ip,ntr)))
3777         $           /clsigma(indmax(cltrx(ip,ntr)))
3778                call applypfa(PFA,cltrx(ip,ntr),angx,corr,res)
3779                xpu(ip,ntr)      = corr
3780    
3781                multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
3782         $                         +10000*mult(cltry(ip,ntr))
3783                seedy(ip,ntr)    = clsignal(indmax(cltry(ip,ntr)))
3784         $           /clsigma(indmax(cltry(ip,ntr)))
3785                call applypfa(PFA,cltry(ip,ntr),angy,corr,res)
3786                ypu(ip,ntr)      = corr
3787    
3788           elseif(icl.ne.0)then           elseif(icl.ne.0)then
3789              if(mod(VIEW(icl),2).eq.0)cltrx(ip,ntr)=icl  
3790              if(mod(VIEW(icl),2).eq.1)cltry(ip,ntr)=icl              cl_used(icl) = 1    !tag used clusters          
3791  c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)  
3792                if(mod(VIEW(icl),2).eq.0)then
3793                   cltrx(ip,ntr)=icl
3794    
3795                   xbad(ip,ntr) = nbadstrips(4,icl)
3796    
3797                   if(nsatstrips(icl).gt.0)dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3798    
3799                   multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
3800         $                         +10000*mult(cltrx(ip,ntr))
3801                   seedx(ip,ntr)    = clsignal(indmax(cltrx(ip,ntr)))
3802         $           /clsigma(indmax(cltrx(ip,ntr)))
3803                   call applypfa(PFA,cltrx(ip,ntr),angx,corr,res)
3804                   xpu(ip,ntr)      = corr
3805    
3806                elseif(mod(VIEW(icl),2).eq.1)then
3807                   cltry(ip,ntr)=icl
3808    
3809                   ybad(ip,ntr) = nbadstrips(4,icl)
3810    
3811                   if(nsatstrips(icl).gt.0)dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3812    
3813                   multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
3814         $                         +10000*mult(cltry(ip,ntr))
3815                   seedy(ip,ntr)    = clsignal(indmax(cltry(ip,ntr)))
3816         $           /clsigma(indmax(cltry(ip,ntr)))
3817                   call applypfa(PFA,cltry(ip,ntr),angy,corr,res)
3818                   ypu(ip,ntr)      = corr
3819                  
3820                endif
3821    
3822           endif                     endif          
3823    
3824        enddo        enddo
3825    
3826          if(DEBUG.eq.1)then
3827             print*,'> STORING TRACK ',ntr
3828             print*,'clusters: '
3829             do ip=1,6
3830                print*,'> ',ip,' -- ',cltrx(ip,ntr),cltry(ip,ntr)
3831             enddo
3832          endif
3833    
3834    c$$$      print*,(xgood(i),i=1,6)
3835    c$$$      print*,(ygood(i),i=1,6)
3836    c$$$      print*,(ls(i,ntr),i=1,6)
3837    c$$$      print*,(dedx_x(i,ntr),i=1,6)
3838    c$$$      print*,(dedx_y(i,ntr),i=1,6)
3839    c$$$      print*,'-----------------------'
3840    
3841        end        end
3842    
# Line 3280  c            print*,ip,' ',cltrx(ip,ntr) Line 3849  c            print*,ip,' ',cltrx(ip,ntr)
3849  *     -------------------------------------------------------  *     -------------------------------------------------------
3850    
3851        include 'commontracker.f'        include 'commontracker.f'
 c      include 'level1.f'  
3852        include 'calib.f'        include 'calib.f'
3853        include 'level1.f'        include 'level1.f'
3854        include 'common_momanhough.f'        include 'common_momanhough.f'
# Line 3288  c      include 'level1.f' Line 3856  c      include 'level1.f'
3856        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3857    
3858  *     count #cluster per plane not associated to any track  *     count #cluster per plane not associated to any track
 c      good2=1!.true.  
3859        nclsx = 0        nclsx = 0
3860        nclsy = 0        nclsy = 0
3861    
3862        do iv = 1,nviews        do iv = 1,nviews
3863           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)
3864             good2(iv) = good2(iv) + mask_view(iv)*2**8
3865        enddo        enddo
3866    
3867          if(DEBUG.eq.1)then
3868             print*,'> STORING SINGLETS '
3869          endif
3870    
3871        do icl=1,nclstr1        do icl=1,nclstr1
3872    
3873             ip=nplanes-npl(VIEW(icl))+1            
3874            
3875           if(cl_used(icl).eq.0)then !cluster not included in any track           if(cl_used(icl).eq.0)then !cluster not included in any track
             ip=nplanes-npl(VIEW(icl))+1              
3876              if(mod(VIEW(icl),2).eq.0)then !=== X views              if(mod(VIEW(icl),2).eq.0)then !=== X views
3877                 nclsx = nclsx + 1                 nclsx = nclsx + 1
3878                 planex(nclsx) = ip                 planex(nclsx) = ip
3879                 sgnlxs(nclsx) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3880                   if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx)
3881                 clsx(nclsx)   = icl                 clsx(nclsx)   = icl
3882                 do is=1,2                 do is=1,2
3883  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
3884                    call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)  c                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
3885                      call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.,0.,0.)
3886                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2
3887                 enddo                 enddo
3888  c$$$               print*,'nclsx         ',nclsx  c$$$               print*,'nclsx         ',nclsx
# Line 3317  c$$$               print*,'xs(2,nclsx)   Line 3893  c$$$               print*,'xs(2,nclsx)  
3893              else                          !=== Y views              else                          !=== Y views
3894                 nclsy = nclsy + 1                 nclsy = nclsy + 1
3895                 planey(nclsy) = ip                 planey(nclsy) = ip
3896                 sgnlys(nclsy) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3897                   if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy)
3898                 clsy(nclsy)   = icl                 clsy(nclsy)   = icl
3899                 do is=1,2                 do is=1,2
3900  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
3901                    call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)  c                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
3902                      call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.,0.,0.)
3903                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2
3904                 enddo                 enddo
3905  c$$$               print*,'nclsy         ',nclsy  c$$$               print*,'nclsy         ',nclsy
# Line 3331  c$$$               print*,'ys(1,nclsy)   Line 3909  c$$$               print*,'ys(1,nclsy)  
3909  c$$$               print*,'ys(2,nclsy)   ',ys(2,nclsy)  c$$$               print*,'ys(2,nclsy)   ',ys(2,nclsy)
3910              endif              endif
3911           endif           endif
 c      print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip)  
3912    
3913  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO
3914           whichtrack(icl) = cl_used(icl)           whichtrack(icl) = cl_used(icl)
3915    *     --------------------------------------------------
3916    *     per non perdere la testa...
3917    *     whichtrack(icl) e` una variabile del common level1
3918    *     che serve solo per sapere quali cluster sono stati
3919    *     associati ad una traccia, e permettere di salvare
3920    *     solo questi nell'albero di uscita
3921    *     --------------------------------------------------
3922            
3923    
3924    c$$$         print*,' cl ',icl,' --> ',cl_used(icl)
3925    c$$$
3926    c$$$         if( cl_used(icl).ne.0 )then
3927    c$$$            if(
3928    c$$$     $           mod(VIEW(icl),2).eq.0.and.
3929    c$$$     $           cltrx(ip,whichtrack(icl)).ne.icl )
3930    c$$$     $           print*,'**WARNING** cltrx(',ip,',',whichtrack(icl)
3931    c$$$     $           ,')=',cltrx(ip,whichtrack(icl)),'.ne.',icl
3932    c$$$            if(
3933    c$$$     $           mod(VIEW(icl),2).eq.1.and.
3934    c$$$     $           cltry(ip,whichtrack(icl)).ne.icl )
3935    c$$$     $           print*,'**WARNING** cltry(',ip,',',whichtrack(icl)
3936    c$$$     $           ,')=',cltry(ip,whichtrack(icl)),'.ne.',icl
3937    c$$$         endif
3938            
3939    
3940        enddo        enddo
3941        end        end
Legend:
Removed from v.1.16  
changed lines
  Added in v.1.32
  ViewVC Help
Powered by ViewVC 1.1.23