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