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

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

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

revision 1.16 by pam-fi, Thu Nov 30 17:04:27 2006 UTC revision 1.28 by pam-fi, Mon Aug 20 16:07:16 2007 UTC
# Line 20  Line 20 
20        include 'calib.f'        include 'calib.f'
21        include 'level2.f'        include 'level2.f'
22    
23  c      include 'momanhough_init.f'  
24    
25    c      print*,'======================================================'
26    c$$$      do ic=1,NCLSTR1
27    c$$$         if(.false.
28    c$$$     $        .or.nsatstrips(ic).gt.0
29    c$$$c     $        .or.nbadstrips(0,ic).gt.0
30    c$$$c     $        .or.nbadstrips(4,ic).gt.0
31    c$$$c     $        .or.nbadstrips(3,ic).gt.0
32    c$$$     $        .or..false.)then
33    c$$$            print*,'--- cl-',ic,' ------------------------'
34    c$$$            istart = INDSTART(IC)
35    c$$$            istop  = TOTCLLENGTH
36    c$$$            if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
37    c$$$            print*,'ADC   ',(CLADC(i),i=istart,istop)
38    c$$$            print*,'s/n   ',(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
39    c$$$            print*,'sgnl  ',(CLSIGNAL(i),i=istart,istop)
40    c$$$            print*,'strip ',(i-INDMAX(ic),i=istart,istop)
41    c$$$            print*,'view ',VIEW(ic)
42    c$$$            print*,'maxs ',MAXS(ic)
43    c$$$            print*,'COG4 ',cog(4,ic)
44    c$$$            ff = fbad_cog(4,ic)
45    c$$$            print*,'fbad ',ff
46    c$$$            print*,(CLBAD(i),i=istart,istop)
47    c$$$            bb=nbadstrips(0,ic)
48    c$$$            print*,'#BAD (tot)',bb
49    c$$$            bb=nbadstrips(4,ic)
50    c$$$            print*,'#BAD (4)',bb
51    c$$$            bb=nbadstrips(3,ic)
52    c$$$            print*,'#BAD (3)',bb
53    c$$$            ss=nsatstrips(ic)
54    c$$$            print*,'#saturated ',ss
55    c$$$         endif
56    c$$$      enddo
57                
58  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
59  *     STEP 1  *     STEP 1
# Line 41  c      include 'momanhough_init.f' Line 74  c      include 'momanhough_init.f'
74  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
75  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
76    
77  c      iflag=0  
78        call cl_to_couples(iflag)        call cl_to_couples(iflag)
79        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
80           goto 880               !go to next event           goto 880               !go to next event
# Line 75  c      iflag=0 Line 108  c      iflag=0
108  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
109  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
110    
111  c      iflag=0  
112        call cp_to_doubtrip(iflag)        call cp_to_doubtrip(iflag)
113        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
114           goto 880               !go to next event                       goto 880               !go to next event            
# Line 182  c      iflag=0 Line 215  c      iflag=0
215  c      include 'momanhough_init.f'  c      include 'momanhough_init.f'
216                
217        logical FIMAGE            !        logical FIMAGE            !
218          real trackimage(NTRACKSMAX)
219        real*8 AL_GUESS(5)        real*8 AL_GUESS(5)
220    
221  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
# Line 239  c$$$         if(ibest.eq.0)goto 880 !>> Line 273  c$$$         if(ibest.eq.0)goto 880 !>>
273  *     2nd) increasing chi**2  *     2nd) increasing chi**2
274  *     -------------------------------------------------------  *     -------------------------------------------------------
275           rchi2best=1000000000.           rchi2best=1000000000.
276           ndofbest=0             !(1)           ndofbest=0            
277           do i=1,ntracks           do i=1,ntracks
278             ndof=0               !(1)             ndof=0              
279             do ii=1,nplanes      !(1)             do ii=1,nplanes    
280               ndof=ndof          !(1)               ndof=ndof        
281       $            +int(xgood_store(ii,i)) !(1)       $            +int(xgood_store(ii,i))
282       $            +int(ygood_store(ii,i)) !(1)       $            +int(ygood_store(ii,i))
283             enddo                !(1)             enddo              
284             if(ndof.gt.ndofbest)then !(1)             if(ndof.gt.ndofbest)then
285               ibest=i               ibest=i
286               rchi2best=RCHI2_STORE(i)               rchi2best=RCHI2_STORE(i)
287               ndofbest=ndof      !(1)               ndofbest=ndof    
288             elseif(ndof.eq.ndofbest)then !(1)             elseif(ndof.eq.ndofbest)then
289               if(RCHI2_STORE(i).lt.rchi2best.and.               if(RCHI2_STORE(i).lt.rchi2best.and.
290       $            RCHI2_STORE(i).gt.0)then       $            RCHI2_STORE(i).gt.0)then
291                 ibest=i                 ibest=i
292                 rchi2best=RCHI2_STORE(i)                 rchi2best=RCHI2_STORE(i)
293                 ndofbest=ndof    !(1)                 ndofbest=ndof  
294               endif              !(1)               endif            
295             endif             endif
296           enddo           enddo
297    
# Line 298  c$$$         enddo Line 332  c$$$         enddo
332              iimage=0              iimage=0
333           endif           endif
334           if(icand.eq.0)then           if(icand.eq.0)then
335              print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand              if(VERBOSE.EQ.1)then
336       $           ,ibest,iimage                 print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand
337         $              ,ibest,iimage
338                endif
339              return              return
340           endif           endif
341    
# Line 314  c$$$         enddo Line 350  c$$$         enddo
350           do i=1,5           do i=1,5
351              AL_GUESS(i)=AL(i)              AL_GUESS(i)=AL(i)
352           enddo           enddo
353    c         print*,'## guess: ',al
354    
355           do i=1,5           do i=1,5
356              AL(i)=dble(AL_STORE(i,icand))                          AL(i)=dble(AL_STORE(i,icand))            
# Line 324  c$$$         enddo Line 361  c$$$         enddo
361           jstep=0                !# minimization steps           jstep=0                !# minimization steps
362    
363           iprint=0           iprint=0
364  c         if(DEBUG)iprint=1  c         if(DEBUG.EQ.1)iprint=1
365           if(VERBOSE)iprint=1           if(VERBOSE.EQ.1)iprint=1
366           if(DEBUG)iprint=2           if(DEBUG.EQ.1)iprint=2
367           call mini2(jstep,ifail,iprint)           call mini2(jstep,ifail,iprint)
368           if(ifail.ne.0) then           if(ifail.ne.0) then
369              if(VERBOSE)then              if(VERBOSE.EQ.1)then
370                 print *,                 print *,
371       $              '*** MINIMIZATION FAILURE *** (after refinement) '       $              '*** MINIMIZATION FAILURE *** (after refinement) '
372       $              ,iev       $              ,iev
# Line 344  c$$$               print*,'------------- Line 381  c$$$               print*,'-------------
381  c            chi2=-chi2  c            chi2=-chi2
382           endif           endif
383                    
384           if(DEBUG)then           if(DEBUG.EQ.1)then
385              print*,'----------------------------- improved track coord'              print*,'----------------------------- improved track coord'
386  22222       format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)  22222       format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
387              do ip=1,6              do ip=1,6
# Line 355  c            chi2=-chi2 Line 392  c            chi2=-chi2
392           endif           endif
393    
394  c         rchi2=chi2/dble(ndof)  c         rchi2=chi2/dble(ndof)
395           if(DEBUG)then           if(DEBUG.EQ.1)then
396              print*,' '              print*,' '
397              print*,'****** SELECTED TRACK *************'              print*,'****** SELECTED TRACK *************'
398              print*,'#         R. chi2        RIG'              print*,'#         R. chi2        RIG'
# Line 371  c         rchi2=chi2/dble(ndof) Line 408  c         rchi2=chi2/dble(ndof)
408  *     ------------- search if the track has an IMAGE -------------  *     ------------- search if the track has an IMAGE -------------
409  *     ------------- (also this is stored )           -------------  *     ------------- (also this is stored )           -------------
410           if(FIMAGE)goto 122     !>>> jump! (this is already an image)           if(FIMAGE)goto 122     !>>> jump! (this is already an image)
411  *     now search for track-image, by comparing couples IDs  
412    *     -----------------------------------------------------
413    *     first check if the track is ambiguous
414    *     -----------------------------------------------------
415    *     (modified on august 2007 by ElenaV)
416             is1=0
417             do ip=1,NPLANES
418                if(SENSOR_STORE(ip,icand).ne.0)then
419                   is1=SENSOR_STORE(ip,icand)
420                   if(ip.eq.6)is1=3-is1 !last plane inverted
421                endif
422             enddo
423             if(is1.eq.0)then
424                if(WARNING.EQ.1)print*,'** WARNING ** sensor=0'
425                goto 122            !jump
426             endif
427    c         print*,'is1 ',is1
428             do ip=1,NPLANES
429                is2 = SENSOR_STORE(ip,icand) !sensor
430    c            print*,'is2 ',is2,' ip ',ip
431                if(ip.eq.6.and.is2.ne.0)is2=3-is2 !last plane inverted
432                if(
433         $           (is1.ne.is2.and.is2.ne.0)
434         $           )goto 122      !jump (this track cannot have an image)
435             enddo
436             if(DEBUG.eq.1)print*,' >>> ambiguous track! '
437    *     now search for track-image among track candidates
438    c$$$         do i=1,ntracks
439    c$$$            iimage=i
440    c$$$            do ip=1,nplanes
441    c$$$               if(     CP_STORE(nplanes-ip+1,icand).ne.
442    c$$$     $              -1*CP_STORE(nplanes-ip+1,i).and.
443    c$$$     $              CP_STORE(nplanes-ip+1,i).ne.0.and.
444    c$$$     $              CP_STORE(nplanes-ip+1,icand).ne.0 )iimage=0
445    c$$$               print*,' track ',i,' CP ',CP_STORE(nplanes-ip+1,i)
446    c$$$     $              ,CP_STORE(nplanes-ip+1,icand),iimage
447    c$$$            enddo
448    c$$$            if(  iimage.ne.0.and.
449    c$$$c     $           RCHI2_STORE(i).le.CHI2MAX.and.
450    c$$$c     $           RCHI2_STORE(i).gt.0.and.
451    c$$$     $           .true.)then
452    c$$$               if(DEBUG.EQ.1)print*,'Track candidate ',iimage
453    c$$$     $              ,' >>> TRACK IMAGE >>> of'
454    c$$$     $              ,ibest
455    c$$$               goto 122         !image track found
456    c$$$            endif
457    c$$$         enddo
458    *     ---------------------------------------------------------------
459    *     take the candidate with the greatest number of matching couples
460    *     if more than one satisfies the requirement,
461    *     choose the one with more points and lower chi2
462    *     ---------------------------------------------------------------
463    *     count the number of matching couples
464             ncpmax = 0
465             iimage   = 0           !id of candidate with better matching
466           do i=1,ntracks           do i=1,ntracks
467              iimage=i              ncp=0
468              do ip=1,nplanes              do ip=1,nplanes
469                 if(     CP_STORE(nplanes-ip+1,icand).ne.                 if(CP_STORE(nplanes-ip+1,icand).ne.0)then
470       $              -1*CP_STORE(nplanes-ip+1,i) )iimage=0                    if(
471         $                 CP_STORE(nplanes-ip+1,i).ne.0
472         $                 .and.
473         $                 CP_STORE(nplanes-ip+1,icand).eq.
474         $                 -1*CP_STORE(nplanes-ip+1,i)
475         $                 )then
476                         ncp=ncp+1  !ok
477                      elseif(
478         $                    CP_STORE(nplanes-ip+1,i).ne.0
479         $                    .and.
480         $                    CP_STORE(nplanes-ip+1,icand).ne.
481         $                    -1*CP_STORE(nplanes-ip+1,i)
482         $                    )then
483                         ncp = 0
484                         goto 117   !it is not an image candidate
485                      else
486                        
487                      endif
488                   endif
489    c$$$               print*,' track ',i,' CP ',CP_STORE(nplanes-ip+1,i)
490    c$$$     $              ,CP_STORE(nplanes-ip+1,icand),ncp
491              enddo              enddo
492              if(  iimage.ne.0.and.   117        continue
493  c     $           RCHI2_STORE(i).le.CHI2MAX.and.              trackimage(i)=ncp   !number of matching couples
494  c     $           RCHI2_STORE(i).gt.0.and.              if(ncp>ncpmax)then
495       $           .true.)then                 ncpmax=ncp
496                 if(DEBUG)print*,'Track candidate ',iimage                 iimage=i
      $              ,' >>> TRACK IMAGE >>> of'  
      $              ,ibest  
                goto 122         !image track found  
497              endif              endif
498           enddo           enddo
499    *     check if there are more than one image candidates
500             ngood=0
501             do i=1,ntracks
502                if( ncpmax.ne.0.and.trackimage(i).eq.ncpmax )ngood=ngood+1
503             enddo
504    *     if there are, choose the best one
505             if(ngood.gt.1)then
506    *     -------------------------------------------------------
507    *     order track-candidates according to:
508    *     1st) decreasing n.points
509    *     2nd) increasing chi**2
510    *     -------------------------------------------------------
511                rchi2best=1000000000.
512                ndofbest=0            
513                do i=1,ntracks
514                   if( trackimage(i).eq.ncpmax )then
515                      ndof=0              
516                      do ii=1,nplanes    
517                         ndof=ndof        
518         $                    +int(xgood_store(ii,i))
519         $                    +int(ygood_store(ii,i))
520                      enddo              
521                      if(ndof.gt.ndofbest)then
522                         iimage=i
523                         rchi2best=RCHI2_STORE(i)
524                         ndofbest=ndof    
525                      elseif(ndof.eq.ndofbest)then
526                         if(RCHI2_STORE(i).lt.rchi2best.and.
527         $                    RCHI2_STORE(i).gt.0)then
528                            iimage=i
529                            rchi2best=RCHI2_STORE(i)
530                            ndofbest=ndof  
531                         endif            
532                      endif
533                   endif
534                enddo
535                
536             endif
537    
538             if(DEBUG.EQ.1)print*,'Track candidate ',iimage
539         $        ,' >>> TRACK IMAGE >>> of'
540         $        ,ibest
541    
542   122     continue   122     continue
543    
544    
545  *     --- and store the results --------------------------------  *     --- and store the results --------------------------------
546           ntrk = ntrk + 1                   !counter of found tracks           ntrk = ntrk + 1                   !counter of found tracks
547           if(.not.FIMAGE           if(.not.FIMAGE
# Line 402  c         print*,'++++++++++ iimage,fima Line 554  c         print*,'++++++++++ iimage,fima
554  c     $        ,iimage,fimage,ntrk,image(ntrk)  c     $        ,iimage,fimage,ntrk,image(ntrk)
555    
556           if(ntrk.eq.NTRKMAX)then           if(ntrk.eq.NTRKMAX)then
557              if(verbose)              if(verbose.eq.1)
558       $           print*,       $           print*,
559       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
560       $           'tracks exceeds vector dimension '       $           'tracks exceeds vector dimension '
# Line 438  cc            good2=.false. Line 590  cc            good2=.false.
590       $        rchi2best.le.CHI2MAX.and.       $        rchi2best.le.CHI2MAX.and.
591  c     $        rchi2best.lt.15..and.  c     $        rchi2best.lt.15..and.
592       $        .true.)then       $        .true.)then
593              if(DEBUG)then              if(DEBUG.EQ.1)then
594                 print*,'***** NEW SEARCH ****'                 print*,'***** NEW SEARCH ****'
595              endif              endif
596              goto 11111          !try new search              goto 11111          !try new search
# Line 476  c     $        rchi2best.lt.15..and. Line 628  c     $        rchi2best.lt.15..and.
628  *     PFAy   - Position Finding Algorithm in y (COG2,ETA2,...)  *     PFAy   - Position Finding Algorithm in y (COG2,ETA2,...)
629  *     angx   - Projected angle in x  *     angx   - Projected angle in x
630  *     angy   - Projected angle in y  *     angy   - Projected angle in y
631    *     bfx    - x-component of magnetci field
632    *     bfy    - y-component of magnetic field
633  *  *
634  *     --------- COUPLES -------------------------------------------------------  *     --------- COUPLES -------------------------------------------------------
635  *     The couple defines a point in the space.  *     The couple defines a point in the space.
# Line 514  c     $        rchi2best.lt.15..and. Line 668  c     $        rchi2best.lt.15..and.
668  *  *
669  *  *
670    
671        subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,angx,angy)        subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
672    
 c*****************************************************  
 c     07/10/2005 modified by elena vannuccini --> (1)  
 c     01/02/2006 modified by elena vannuccini --> (2)  
 c     02/02/2006 modified by Elena Vannuccini --> (3)  
 c                (implemented new p.f.a.)  
 c     03/02/2006 modified by Elena Vannuccini --> (4)  
 c                (implemented variable resolution)  
 c*****************************************************  
673                
674        include 'commontracker.f'        include 'commontracker.f'
675        include 'level1.f'        include 'level1.f'
676        include 'calib.f'        include 'calib.f'
 c      include 'level1.f'  
677        include 'common_align.f'        include 'common_align.f'
678        include 'common_mech.f'        include 'common_mech.f'
679        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
 c      include 'common_resxy.f'  
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
680    
681        integer icx,icy           !X-Y cluster ID        integer icx,icy           !X-Y cluster ID
682        integer sensor        integer sensor
683        integer viewx,viewy        integer viewx,viewy
684        character*4 PFAx,PFAy     !PFA to be used        character*4 PFAx,PFAy     !PFA to be used
685        real angx,angy            !X-Y angle        real ax,ay                !X-Y geometric angle
686          real angx,angy            !X-Y effective angle
687          real bfx,bfy              !X-Y b-field components
688    
689        real stripx,stripy        real stripx,stripy
690    
691        double precision xrt,yrt,zrt        double precision xrt,yrt,zrt
692        double precision xrt_A,yrt_A,zrt_A        double precision xrt_A,yrt_A,zrt_A
693        double precision xrt_B,yrt_B,zrt_B        double precision xrt_B,yrt_B,zrt_B
 c      double precision xi,yi,zi  
 c      double precision xi_A,yi_A,zi_A  
 c      double precision xi_B,yi_B,zi_B  
694                
695    
696        parameter (ndivx=30)        parameter (ndivx=30)
697    
698    
699    c$$$      print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy
700                
701        resxPAM = 0        resxPAM = 0
702        resyPAM = 0        resyPAM = 0
# Line 567  c      double precision xi_B,yi_B,zi_B Line 710  c      double precision xi_B,yi_B,zi_B
710        xPAM_B = 0.        xPAM_B = 0.
711        yPAM_B = 0.        yPAM_B = 0.
712        zPAM_B = 0.        zPAM_B = 0.
713    c      print*,'## xyz_PAM: ',icx,icy,sensor,PFAx,PFAy,angx,angy
714    
715          if(sensor.lt.1.or.sensor.gt.2)then
716             print*,'xyz_PAM   ***ERROR*** wrong input '
717             print*,'sensor ',sensor
718             icx=0
719             icy=0
720          endif
721    
722  *     -----------------  *     -----------------
723  *     CLUSTER X  *     CLUSTER X
724  *     -----------------  *     -----------------      
   
725        if(icx.ne.0)then        if(icx.ne.0)then
726           viewx = VIEW(icx)  
727           nldx = nld(MAXS(icx),VIEW(icx))           viewx   = VIEW(icx)
728           nplx = npl(VIEW(icx))           nldx    = nld(MAXS(icx),VIEW(icx))
729           resxPAM = RESXAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!           nplx    = npl(VIEW(icx))
730                     resxPAM = RESXAV
731           stripx = float(MAXS(icx))           stripx  = float(MAXS(icx))
732           if(PFAx.eq.'COG1')then  !(1)  
733              stripx = stripx      !(1)           if(
734              resxPAM = resxPAM    !(1)       $        viewx.lt.1.or.        
735         $        viewx.gt.12.or.
736         $        nldx.lt.1.or.
737         $        nldx.gt.3.or.
738         $        stripx.lt.1.or.
739         $        stripx.gt.3072.or.
740         $        .false.)then
741                print*,'xyz_PAM   ***ERROR*** wrong input '
742                print*,'icx ',icx,'view ',viewx,'nld ',nldx,'strip ',stripx
743                icx = 0
744                goto 10
745             endif
746    
747    *        --------------------------
748    *        magnetic-field corrections
749    *        --------------------------
750             angtemp  = ax
751             bfytemp  = bfy
752    *        /////////////////////////////////
753    *        AAAAHHHHHHHH!!!!!!!!!!!!!!!!!!!!!
754    *        *grvzkkjsdgjhhhgngbn###>:(
755    *        /////////////////////////////////
756    c         if(nplx.eq.6) angtemp = -1. * ax
757    c         if(nplx.eq.6) bfytemp = -1. * bfy
758             if(viewx.eq.12) angtemp = -1. * ax
759             if(viewx.eq.12) bfytemp = -1. * bfy
760             tgtemp   = tan(angtemp*acos(-1.)/180.) + pmuH_h*bfytemp*0.00001
761             angx     = 180.*atan(tgtemp)/acos(-1.)
762             stripx   = stripx - 0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX
763    c$$$         print*,nplx,ax,bfy/10.
764    c$$$         print*,angx,0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX
765    c$$$         print*,'========================'
766    c$$$         if(bfy.ne.0.)print*,viewx,'-x- '
767    c$$$     $        ,bfy,-1*0.5*pmuH_h*bfytemp*0.00001*SiDimZ
768    *        --------------------------
769    
770    c$$$         print*,'--- x-cl ---'
771    c$$$         istart = INDSTART(ICX)
772    c$$$         istop  = TOTCLLENGTH
773    c$$$         if(icx.lt.NCLSTR1)istop=INDSTART(ICX+1)-1
774    c$$$         print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
775    c$$$         print*,(CLSIGNAL(i),i=istart,istop)
776    c$$$         print*,INDMAX(icx)
777    c$$$         print*,cog(4,icx)
778    c$$$         print*,fbad_cog(4,icx)
779            
780    
781             if(PFAx.eq.'COG1')then
782    
783                stripx  = stripx
784                resxPAM = 1e-4*pitchX/sqrt(12.)!!resxPAM
785    
786           elseif(PFAx.eq.'COG2')then           elseif(PFAx.eq.'COG2')then
787              stripx = stripx + cog(2,icx)              
788                stripx  = stripx + cog(2,icx)            
789                resxPAM = risx_cog(abs(angx))!TEMPORANEO              
790              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
791    
792             elseif(PFAx.eq.'COG3')then
793    
794                stripx  = stripx + cog(3,icx)            
795                resxPAM = risx_cog(abs(angx))!TEMPORANEO                      
796                resxPAM = resxPAM*fbad_cog(3,icx)
797    
798             elseif(PFAx.eq.'COG4')then
799    
800                stripx  = stripx + cog(4,icx)            
801                resxPAM = risx_cog(abs(angx))!TEMPORANEO                      
802                resxPAM = resxPAM*fbad_cog(4,icx)
803    
804           elseif(PFAx.eq.'ETA2')then           elseif(PFAx.eq.'ETA2')then
805  c            cog2 = cog(2,icx)  
806  c            etacorr = pfaeta2(cog2,viewx,nldx,angx)                          stripx  = stripx + pfaeta2(icx,angx)          
807  c            stripx = stripx + etacorr              resxPAM = risxeta2(abs(angx))
             stripx = stripx + pfaeta2(icx,angx)            !(3)  
             resxPAM = risx_eta2(angx)                       !   (4)  
             if(DEBUG.and.fbad_cog(2,icx).ne.1)  
      $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)  
808              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
809           elseif(PFAx.eq.'ETA3')then                         !(3)  c$$$            if(DEBUG.EQ.1.and.fbad_cog(2,icx).ne.1)
810              stripx = stripx + pfaeta3(icx,angx)            !(3)  c$$$     $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
811              resxPAM = risx_eta3(angx)                       !   (4)  
812              if(DEBUG.and.fbad_cog(3,icx).ne.1)              !(3)           elseif(PFAx.eq.'ETA3')then                        
813       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)  
814              resxPAM = resxPAM*fbad_cog(3,icx)               !(3)              stripx  = stripx + pfaeta3(icx,angx)          
815           elseif(PFAx.eq.'ETA4')then                         !(3)              resxPAM = risxeta3(abs(angx))                      
816              stripx = stripx + pfaeta4(icx,angx)            !(3)              resxPAM = resxPAM*fbad_cog(3,icx)              
817              resxPAM = risx_eta4(angx)                       !   (4)  c            if(DEBUG.and.fbad_cog(3,icx).ne.1)            
818              if(DEBUG.and.fbad_cog(4,icx).ne.1)              !(3)  c     $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)
819       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)  
820              resxPAM = resxPAM*fbad_cog(4,icx)               !(3)           elseif(PFAx.eq.'ETA4')then                        
821           elseif(PFAx.eq.'ETA')then                          !(3)  
822              stripx = stripx + pfaeta(icx,angx)             !(3)              stripx  = stripx + pfaeta4(icx,angx)            
823              resxPAM = ris_eta(icx,angx)                     !   (4)              resxPAM = risxeta4(abs(angx))                      
824              if(DEBUG.and.fbad_cog(2,icx).ne.1)              !(3)              resxPAM = resxPAM*fbad_cog(4,icx)              
825       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)  c            if(DEBUG.and.fbad_cog(4,icx).ne.1)              
826  c            resxPAM = resxPAM*fbad_cog(2,icx)               !(3)TEMPORANEO  c     $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)
827              resxPAM = resxPAM*fbad_eta(icx,angx)             !(3)(4)  
828           elseif(PFAx.eq.'COG')then           !(2)           elseif(PFAx.eq.'ETA')then  
829              stripx = stripx + cog(0,icx)     !(2)          
830              resxPAM = risx_cog(angx)                        !   (4)              stripx  = stripx + pfaeta(icx,angx)            
831              resxPAM = resxPAM*fbad_cog(0,icx)!(2)  c            resxPAM = riseta(icx,angx)                    
832                resxPAM = riseta(viewx,angx)                    
833                resxPAM = resxPAM*fbad_eta(icx,angx)            
834    c            if(DEBUG.and.fbad_cog(2,icx).ne.1)              
835    c     $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
836    
837             elseif(PFAx.eq.'ETAL')then  
838    
839                stripx  = stripx + pfaetal(icx,angx)            
840                resxPAM = riseta(viewx,angx)                    
841                resxPAM = resxPAM*fbad_eta(icx,angx)            
842    c            if(DEBUG.and.fbad_cog(2,icx).ne.1)              
843    c     $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
844    
845             elseif(PFAx.eq.'COG')then          
846    
847                stripx  = stripx + cog(0,icx)            
848                resxPAM = risx_cog(abs(angx))                    
849                resxPAM = resxPAM*fbad_cog(0,icx)
850    
851           else           else
852              print*,'*** Non valid p.f.a. (x) --> ',PFAx              if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (x) --> ',PFAx
853           endif           endif
854    
855        endif  
856  c      if(icy.eq.0.and.icx.ne.0)  *     ======================================
857  c     $     print*,PFAx,icx,angx,stripx,resxPAM,'***'  *     temporary patch for saturated clusters
858          *     ======================================
859             if( nsatstrips(icx).gt.0 )then
860                stripx  = stripx + cog(4,icx)            
861                resxPAM = pitchX*1e-4/sqrt(12.)
862                cc=cog(4,icx)
863    c$$$            print*,icx,' *** ',cc
864    c$$$            print*,icx,' *** ',resxPAM
865             endif
866    
867     10   endif
868    
869        
870  *     -----------------  *     -----------------
871  *     CLUSTER Y  *     CLUSTER Y
872  *     -----------------  *     -----------------
873    
874        if(icy.ne.0)then        if(icy.ne.0)then
875    
876           viewy = VIEW(icy)           viewy = VIEW(icy)
877           nldy = nld(MAXS(icy),VIEW(icy))           nldy = nld(MAXS(icy),VIEW(icy))
878           nply = npl(VIEW(icy))           nply = npl(VIEW(icy))
879           resyPAM = RESYAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!           resyPAM = RESYAV
880             stripy = float(MAXS(icy))
881    
882             if(
883         $        viewy.lt.1.or.        
884         $        viewy.gt.12.or.
885         $        nldy.lt.1.or.
886         $        nldy.gt.3.or.
887         $        stripy.lt.1.or.
888         $        stripy.gt.3072.or.
889         $        .false.)then
890                print*,'xyz_PAM   ***ERROR*** wrong input '
891                print*,'icy ',icy,'view ',viewy,'nld ',nldy,'strip ',stripy
892                icy = 0
893                goto 20
894             endif
895    
896           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
897              print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '              if(DEBUG.EQ.1) then
898       $           ,icx,icy                 print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '
899         $              ,icx,icy
900                endif
901              goto 100              goto 100
902           endif           endif
903    *        --------------------------
904    *        magnetic-field corrections
905    *        --------------------------
906             tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
907             angy    = 180.*atan(tgtemp)/acos(-1.)
908             stripy = stripy + 0.5*pmuH_e*bfx*0.00001*SiDimZ/pitchY
909    c$$$         if(bfx.ne.0.)print*,viewy,'-y- '
910    c$$$     $        ,bfx,0.5*pmuH_e*bfx*0.00001*SiDimZ
911    *        --------------------------
912                    
913           stripy = float(MAXS(icy))  c$$$         print*,'--- y-cl ---'
914           if(PFAy.eq.'COG1')then !(1)  c$$$         istart = INDSTART(ICY)
915              stripy = stripy     !(1)  c$$$         istop  = TOTCLLENGTH
916              resyPAM = resyPAM   !(1)  c$$$         if(icy.lt.NCLSTR1)istop=INDSTART(ICY+1)-1
917    c$$$         print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
918    c$$$         print*,(CLSIGNAL(i),i=istart,istop)
919    c$$$         print*,INDMAX(icy)
920    c$$$         print*,cog(4,icy)
921    c$$$         print*,fbad_cog(4,icy)
922    
923             if(PFAy.eq.'COG1')then
924    
925                stripy  = stripy    
926                resyPAM = 1e-4*pitchY/sqrt(12.)!resyPAM  
927    
928           elseif(PFAy.eq.'COG2')then           elseif(PFAy.eq.'COG2')then
929              stripy = stripy + cog(2,icy)  
930                stripy  = stripy + cog(2,icy)
931                resyPAM = risy_cog(abs(angy))!TEMPORANEO
932              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
933    
934             elseif(PFAy.eq.'COG3')then
935    
936                stripy  = stripy + cog(3,icy)
937                resyPAM = risy_cog(abs(angy))!TEMPORANEO
938                resyPAM = resyPAM*fbad_cog(3,icy)
939    
940             elseif(PFAy.eq.'COG4')then
941    
942                stripy  = stripy + cog(4,icy)
943                resyPAM = risy_cog(abs(angy))!TEMPORANEO
944                resyPAM = resyPAM*fbad_cog(4,icy)
945    
946           elseif(PFAy.eq.'ETA2')then           elseif(PFAy.eq.'ETA2')then
947  c            cog2 = cog(2,icy)  
948  c            etacorr = pfaeta2(cog2,viewy,nldy,angy)              stripy  = stripy + pfaeta2(icy,angy)          
949  c            stripy = stripy + etacorr              resyPAM = risyeta2(abs(angy))              
             stripy = stripy + pfaeta2(icy,angy)            !(3)  
             resyPAM = risy_eta2(angy)                       !   (4)  
950              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
951              if(DEBUG.and.fbad_cog(2,icy).ne.1)  c            if(DEBUG.and.fbad_cog(2,icy).ne.1)
952       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)  c     $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
953           elseif(PFAy.eq.'ETA3')then                         !(3)  
954              stripy = stripy + pfaeta3(icy,angy)            !(3)           elseif(PFAy.eq.'ETA3')then                      
955              resyPAM = resyPAM*fbad_cog(3,icy)               !(3)  
956              if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)              stripy  = stripy + pfaeta3(icy,angy)
957       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)              resyPAM = resyPAM*fbad_cog(3,icy)  
958           elseif(PFAy.eq.'ETA4')then                         !(3)  c            if(DEBUG.and.fbad_cog(3,icy).ne.1)
959              stripy = stripy + pfaeta4(icy,angy)            !(3)  c     $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)
960              resyPAM = resyPAM*fbad_cog(4,icy)               !(3)  
961              if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)           elseif(PFAy.eq.'ETA4')then  
962       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)  
963           elseif(PFAy.eq.'ETA')then                          !(3)              stripy  = stripy + pfaeta4(icy,angy)
964              stripy = stripy + pfaeta(icy,angy)             !(3)              resyPAM = resyPAM*fbad_cog(4,icy)
965              resyPAM = ris_eta(icy,angy)                     !   (4)  c            if(DEBUG.and.fbad_cog(4,icy).ne.1)
966  c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO  c     $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)
967              resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)  
968              if(DEBUG.and.fbad_cog(2,icy).ne.1)              !(3)           elseif(PFAy.eq.'ETA')then
969       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)!(3)  
970                stripy  = stripy + pfaeta(icy,angy)
971    c            resyPAM = riseta(icy,angy)  
972                resyPAM = riseta(viewy,angy)  
973                resyPAM = resyPAM*fbad_eta(icy,angy)
974    c            if(DEBUG.and.fbad_cog(2,icy).ne.1)
975    c     $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
976    
977             elseif(PFAy.eq.'ETAL')then
978    
979                stripy  = stripy + pfaetal(icy,angy)
980                resyPAM = riseta(viewy,angy)  
981                resyPAM = resyPAM*fbad_eta(icy,angy)
982    c            if(DEBUG.and.fbad_cog(2,icy).ne.1)
983    c     $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
984    
985           elseif(PFAy.eq.'COG')then           elseif(PFAy.eq.'COG')then
986              stripy = stripy + cog(0,icy)              
987              resyPAM = risy_cog(angy)                        !   (4)              stripy  = stripy + cog(0,icy)            
988  c            resyPAM = ris_eta(icy,angy)                    !   (4)              resyPAM = risy_cog(abs(angy))
989              resyPAM = resyPAM*fbad_cog(0,icy)              resyPAM = resyPAM*fbad_cog(0,icy)
990    
991           else           else
992              print*,'*** Non valid p.f.a. (x) --> ',PFAx              if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (x) --> ',PFAx
993           endif           endif
994    
       endif  
995    
996          *     ======================================
997    *     temporary patch for saturated clusters
998    *     ======================================
999             if( nsatstrips(icy).gt.0 )then
1000                stripy  = stripy + cog(4,icy)            
1001                resyPAM = pitchY*1e-4/sqrt(12.)
1002                cc=cog(4,icy)
1003    c$$$            print*,icy,' *** ',cc
1004    c$$$            print*,icy,' *** ',resyPAM
1005             endif
1006    
1007    
1008     20   endif
1009    
1010    c$$$      print*,'## stripx,stripy ',stripx,stripy
1011    
1012  c===========================================================  c===========================================================
1013  C     COUPLE  C     COUPLE
1014  C===========================================================  C===========================================================
# Line 697  c     (xi,yi,zi) = mechanical coordinate Line 1019  c     (xi,yi,zi) = mechanical coordinate
1019  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1020           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
1021       $        .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...
1022              print*,'xyz_PAM (couple):',              if(DEBUG.EQ.1) then
1023       $          ' WARNING: false X strip: strip ',stripx                 print*,'xyz_PAM (couple):',
1024         $              ' WARNING: false X strip: strip ',stripx
1025                endif
1026           endif           endif
1027           xi = acoordsi(stripx,viewx)           xi = acoordsi(stripx,viewx)
1028           yi = acoordsi(stripy,viewy)           yi = acoordsi(stripy,viewy)
# Line 790  c            print*,'X-singlet ',icx,npl Line 1114  c            print*,'X-singlet ',icx,npl
1114  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...
1115              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
1116       $           .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...
1117                 print*,'xyz_PAM (X-singlet):',                 if(DEBUG.EQ.1) then
1118       $             ' WARNING: false X strip: strip ',stripx                    print*,'xyz_PAM (X-singlet):',
1119         $                 ' WARNING: false X strip: strip ',stripx
1120                   endif
1121              endif              endif
1122              xi   = acoordsi(stripx,viewx)              xi   = acoordsi(stripx,viewx)
1123    
# Line 813  c            print*,'X-cl ',icx,stripx,' Line 1139  c            print*,'X-cl ',icx,stripx,'
1139  c            print*,yi_A,' <--> ',yi_B  c            print*,yi_A,' <--> ',yi_B
1140    
1141           else           else
1142                if(DEBUG.EQ.1) then
1143              print *,'routine xyz_PAM ---> not properly used !!!'                 print *,'routine xyz_PAM ---> not properly used !!!'
1144              print *,'icx = ',icx                 print *,'icx = ',icx
1145              print *,'icy = ',icy                 print *,'icy = ',icy
1146                endif
1147              goto 100              goto 100
1148                            
1149           endif           endif
# Line 881  c--------------------------------------- Line 1208  c---------------------------------------
1208  c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'  c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'
1209    
1210        else        else
1211                       if(DEBUG.EQ.1) then
1212           print *,'routine xyz_PAM ---> not properly used !!!'              print *,'routine xyz_PAM ---> not properly used !!!'
1213           print *,'icx = ',icx              print *,'icx = ',icx
1214           print *,'icy = ',icy              print *,'icy = ',icy
1215                         endif
1216        endif        endif
1217                    
1218    
1219    c      print*,'## xPAM,yPAM,zPAM       ',xPAM,yPAM,zPAM
1220    c      print*,'## xPAM_A,yPAM_A,zPAM_A ',xPAM_A,yPAM_A,zPAM_A
1221    c      print*,'## xPAM_B,yPAM_B,zPAM_B ',xPAM_B,yPAM_B,zPAM_B
1222    
1223   100  continue   100  continue
1224        end        end
1225    
1226    ************************************************************************
1227    *     Call xyz_PAM subroutine with default PFA and fill the mini2 common.
1228    *     (done to be called from c/c++)
1229    ************************************************************************
1230    
1231          subroutine xyzpam(ip,icx,icy,lad,sensor,ax,ay,bfx,bfy)
1232    
1233          include 'commontracker.f'
1234          include 'level1.f'
1235          include 'common_mini_2.f'
1236          include 'common_xyzPAM.f'
1237          include 'common_mech.f'
1238          include 'calib.f'
1239          
1240    *     flag to chose PFA
1241    c$$$      character*10 PFA
1242    c$$$      common/FINALPFA/PFA
1243    
1244          integer icx,icy           !X-Y cluster ID
1245          integer sensor
1246          character*4 PFAx,PFAy     !PFA to be used
1247          real ax,ay                !X-Y geometric angle
1248          real bfx,bfy              !X-Y b-field components
1249    
1250          ipx=0
1251          ipy=0      
1252          
1253    c$$$      PFAx = 'COG4'!PFA
1254    c$$$      PFAy = 'COG4'!PFA
1255    
1256          if(icx.gt.nclstr1.or.icy.gt.nclstr1)then
1257                print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1258         $           ,' does not exists (nclstr1=',nclstr1,')'
1259                icx = -1*icx
1260                icy = -1*icy
1261                return
1262            
1263          endif
1264          
1265          call idtoc(pfaid,PFAx)
1266          call idtoc(pfaid,PFAy)
1267    
1268    c$$$      call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1269    
1270    c$$$      print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy
1271          
1272          if(icx.ne.0.and.icy.ne.0)then
1273    
1274             ipx=npl(VIEW(icx))
1275             ipy=npl(VIEW(icy))
1276    c$$$         if( (nplanes-ipx+1).ne.ip.or.(nplanes-ipy+1).ne.ip )
1277    c$$$     $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1278    c$$$     $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1279            
1280             if( (nplanes-ipx+1).ne.ip )then
1281                print*,'xyzpam: ***WARNING*** cluster ',icx
1282         $           ,' does not belong to plane: ',ip
1283                icx = -1*icx
1284                return
1285             endif
1286             if( (nplanes-ipy+1).ne.ip )then
1287                print*,'xyzpam: ***WARNING*** cluster ',icy
1288         $           ,' does not belong to plane: ',ip
1289                icy = -1*icy
1290                return
1291             endif
1292    
1293             call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1294    
1295             xgood(ip) = 1.
1296             ygood(ip) = 1.
1297             resx(ip)  = resxPAM
1298             resy(ip)  = resyPAM
1299    
1300             xm(ip) = xPAM
1301             ym(ip) = yPAM
1302             zm(ip) = zPAM
1303             xm_A(ip) = 0.
1304             ym_A(ip) = 0.
1305             xm_B(ip) = 0.
1306             ym_B(ip) = 0.
1307    
1308    c         zv(ip) = zPAM
1309    
1310          elseif(icx.eq.0.and.icy.ne.0)then
1311    
1312             ipy=npl(VIEW(icy))
1313    c$$$         if((nplanes-ipy+1).ne.ip)
1314    c$$$     $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1315    c$$$     $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1316             if( (nplanes-ipy+1).ne.ip )then
1317                print*,'xyzpam: ***WARNING*** cluster ',icy
1318         $           ,' does not belong to plane: ',ip
1319                icy = -1*icy
1320                return
1321             endif
1322    
1323             call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1324            
1325             xgood(ip) = 0.
1326             ygood(ip) = 1.
1327             resx(ip)  = 1000.
1328             resy(ip)  = resyPAM
1329    
1330             xm(ip) = -100.
1331             ym(ip) = -100.
1332             zm(ip) = (zPAM_A+zPAM_B)/2.
1333             xm_A(ip) = xPAM_A
1334             ym_A(ip) = yPAM_A
1335             xm_B(ip) = xPAM_B
1336             ym_B(ip) = yPAM_B
1337    
1338    c         zv(ip) = (zPAM_A+zPAM_B)/2.
1339            
1340          elseif(icx.ne.0.and.icy.eq.0)then
1341    
1342             ipx=npl(VIEW(icx))
1343    c$$$         if((nplanes-ipx+1).ne.ip)
1344    c$$$     $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1345    c$$$     $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1346    
1347             if( (nplanes-ipx+1).ne.ip )then
1348                print*,'xyzpam: ***WARNING*** cluster ',icx
1349         $           ,' does not belong to plane: ',ip
1350                icx = -1*icx
1351                return
1352             endif
1353    
1354             call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1355          
1356             xgood(ip) = 1.
1357             ygood(ip) = 0.
1358             resx(ip)  = resxPAM
1359             resy(ip)  = 1000.
1360    
1361             xm(ip) = -100.
1362             ym(ip) = -100.
1363             zm(ip) = (zPAM_A+zPAM_B)/2.
1364             xm_A(ip) = xPAM_A
1365             ym_A(ip) = yPAM_A
1366             xm_B(ip) = xPAM_B
1367             ym_B(ip) = yPAM_B
1368            
1369    c         zv(ip) = (zPAM_A+zPAM_B)/2.
1370    
1371          else
1372    
1373             il = 2
1374             if(lad.ne.0)il=lad
1375             is = 1
1376             if(sensor.ne.0)is=sensor
1377    c         print*,nplanes-ip+1,il,is
1378    
1379             xgood(ip) = 0.
1380             ygood(ip) = 0.
1381             resx(ip)  = 1000.
1382             resy(ip)  = 1000.
1383    
1384             xm(ip) = -100.
1385             ym(ip) = -100.          
1386             zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1387             xm_A(ip) = 0.
1388             ym_A(ip) = 0.
1389             xm_B(ip) = 0.
1390             ym_B(ip) = 0.
1391    
1392    c         zv(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1393    
1394          endif
1395    
1396          if(DEBUG.EQ.1)then
1397    c         print*,'----------------------------- track coord'
1398    22222    format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
1399             write(*,22222)ip,zm(ip),xm(ip),ym(ip)
1400         $        ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip)
1401         $        ,xgood(ip),ygood(ip),resx(ip),resy(ip)
1402    c$$$         print*,'-----------------------------'
1403          endif
1404          end
1405    
1406  ********************************************************************************  ********************************************************************************
1407  ********************************************************************************  ********************************************************************************
# Line 966  c         print*,'A-(',xPAM_A,yPAM_A,') Line 1477  c         print*,'A-(',xPAM_A,yPAM_A,')
1477           endif                   endif        
1478    
1479           distance=           distance=
1480       $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2       $       ((xmi-XPP)**2+(ymi-YPP)**2)!QUIQUI
1481    cc     $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2
1482           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1483    
1484  c$$$         print*,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b  c$$$         print*,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b
# Line 991  c$$$         print*,' resolution ',re Line 1503  c$$$         print*,' resolution ',re
1503  *     ----------------------  *     ----------------------
1504                    
1505           distance=           distance=
1506       $        ((xPAM-XPP)/resxPAM)**2       $       ((xPAM-XPP))**2 !QUIQUI
1507       $        +       $       +
1508       $        ((yPAM-YPP)/resyPAM)**2       $       ((yPAM-YPP))**2
1509    c$$$     $        ((xPAM-XPP)/resxPAM)**2
1510    c$$$     $        +
1511    c$$$     $        ((yPAM-YPP)/resyPAM)**2
1512           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1513    
1514  c$$$         print*,xPAM,yPAM,zPAM  c$$$         print*,xPAM,yPAM,zPAM
# Line 1002  c$$$         print*,' resolution ',resxP Line 1517  c$$$         print*,' resolution ',resxP
1517                    
1518        else        else
1519                    
1520           print*  c         print*
1521       $        ,' function distance_to ---> wrong usage!!!'  c     $        ,' function distance_to ---> wrong usage!!!'
1522           print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM  c         print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM
1523           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 '
1524       $        ,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
1525        endif          endif  
1526    
1527        distance_to = sngl(distance)        distance_to = sngl(distance)
# Line 1074  c--------------------------------------- Line 1589  c---------------------------------------
1589                 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)                 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
1590       $              .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...
1591  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...
1592                    print*,'whichsensor: ',  c                  print*,'whichsensor: ',
1593       $                ' WARNING: false X strip: strip ',stripx  c     $                ' WARNING: false X strip: strip ',stripx
1594                 endif                 endif
1595                 xi = acoordsi(stripx,viewx)                 xi = acoordsi(stripx,viewx)
1596                 yi = acoordsi(stripy,viewy)                 yi = acoordsi(stripy,viewy)
# Line 1230  c      include 'common_analysis.f' Line 1745  c      include 'common_analysis.f'
1745        is_cp=0        is_cp=0
1746        if(id.lt.0)is_cp=1        if(id.lt.0)is_cp=1
1747        if(id.gt.0)is_cp=2        if(id.gt.0)is_cp=2
1748        if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'  c      if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'
1749    
1750        return        return
1751        end        end
# Line 1329  c      include 'level1.f' Line 1844  c      include 'level1.f'
1844        integer iflag        integer iflag
1845    
1846        integer badseed,badclx,badcly        integer badseed,badclx,badcly
1847        
1848          if(DEBUG.EQ.1)print*,'cl_to_couples:'
1849    
1850  *     init variables  *     init variables
1851        ncp_tot=0        ncp_tot=0
# Line 1359  c      include 'level1.f' Line 1876  c      include 'level1.f'
1876  *     mask views with too many clusters  *     mask views with too many clusters
1877        do iv=1,nviews        do iv=1,nviews
1878           if( ncl_view(iv).gt. nclusterlimit)then           if( ncl_view(iv).gt. nclusterlimit)then
1879              mask_view(iv) = 1  c            mask_view(iv) = 1
1880              if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '              mask_view(iv) = mask_view(iv) + 2**0
1881       $           ,nclusterlimit,' on view ', iv,' --> masked!'              if(DEBUG.EQ.1)
1882         $        print*,' * WARNING * cl_to_couple: n.clusters > '
1883         $        ,nclusterlimit,' on view ', iv,' --> masked!'
1884           endif           endif
1885        enddo        enddo
1886    
# Line 1378  c      include 'level1.f' Line 1897  c      include 'level1.f'
1897  *     ----------------------------------------------------  *     ----------------------------------------------------
1898  *     cut on charge (X VIEW)  *     cut on charge (X VIEW)
1899  *     ----------------------------------------------------  *     ----------------------------------------------------
1900           if(dedx(icx).lt.dedx_x_min)then           if(sgnl(icx).lt.dedx_x_min)then
1901              cl_single(icx)=0              cl_single(icx)=0
1902              goto 10              goto 10
1903           endif           endif
# Line 1428  c     endif Line 1947  c     endif
1947  *     ----------------------------------------------------  *     ----------------------------------------------------
1948  *     cut on charge (Y VIEW)  *     cut on charge (Y VIEW)
1949  *     ----------------------------------------------------  *     ----------------------------------------------------
1950              if(dedx(icy).lt.dedx_y_min)then              if(sgnl(icy).lt.dedx_y_min)then
1951                 cl_single(icy)=0                 cl_single(icy)=0
1952                 goto 20                 goto 20
1953              endif              endif
# Line 1474  c     endif Line 1993  c     endif
1993  *     charge correlation  *     charge correlation
1994  *     (modified to be applied only below saturation... obviously)  *     (modified to be applied only below saturation... obviously)
1995    
1996                 if(  .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx)                 if(  .not.(sgnl(icy).gt.chsaty.and.sgnl(icx).gt.chsatx)
1997       $              .and.       $              .and.
1998       $              .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx)       $              .not.(sgnl(icy).lt.chmipy.and.sgnl(icx).lt.chmipx)
1999       $              .and.       $              .and.
2000       $              (badclx.eq.1.and.badcly.eq.1)       $              (badclx.eq.1.and.badcly.eq.1)
2001       $              .and.       $              .and.
2002       $              .true.)then       $              .true.)then
2003    
2004                    ddd=(dedx(icy)                    ddd=(sgnl(icy)
2005       $                 -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))       $                 -kch(nplx,nldx)*sgnl(icx)-cch(nplx,nldx))
2006                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
2007    
2008  c                  cut = chcut * sch(nplx,nldx)  c                  cut = chcut * sch(nplx,nldx)
2009    
2010                    sss=(kch(nplx,nldx)*dedx(icy)+dedx(icx)                    sss=(kch(nplx,nldx)*sgnl(icy)+sgnl(icx)
2011       $                 -kch(nplx,nldx)*cch(nplx,nldx))       $                 -kch(nplx,nldx)*cch(nplx,nldx))
2012                    sss=sss/sqrt(kch(nplx,nldx)**2+1)                    sss=sss/sqrt(kch(nplx,nldx)**2+1)
2013                    cut = chcut * (16 + sss/50.)                    cut = chcut * (16 + sss/50.)
# Line 1499  c                  cut = chcut * sch(npl Line 2018  c                  cut = chcut * sch(npl
2018                 endif                 endif
2019    
2020                 if(ncp_plane(nplx).gt.ncouplemax)then                 if(ncp_plane(nplx).gt.ncouplemax)then
2021                    if(verbose)print*,                    if(verbose.eq.1)print*,
2022       $                 '** warning ** number of identified '//       $                 '** warning ** number of identified '//
2023       $                 'couples on plane ',nplx,       $                 'couples on plane ',nplx,
2024       $                 'exceeds vector dimention '       $                 'exceeds vector dimention '
2025       $                 ,'( ',ncouplemax,' ) --> masked!'       $                 ,'( ',ncouplemax,' ) --> masked!'
2026                    mask_view(nviewx(nplx)) = 2  c                  mask_view(nviewx(nplx)) = 2
2027                    mask_view(nviewy(nply)) = 2  c                  mask_view(nviewy(nply)) = 2
2028                      mask_view(nviewx(nplx))= mask_view(nviewx(nplx))+ 2**1
2029                      mask_view(nviewy(nply))= mask_view(nviewy(nply))+ 2**1
2030                    goto 10                    goto 10
2031                 endif                 endif
2032                                
# Line 1535  c                  cut = chcut * sch(npl Line 2056  c                  cut = chcut * sch(npl
2056        enddo        enddo
2057                
2058                
2059        if(DEBUG)then        if(DEBUG.EQ.1)then
2060           print*,'clusters  ',nclstr1           print*,'clusters  ',nclstr1
2061           print*,'good    ',(cl_good(i),i=1,nclstr1)           print*,'good    ',(cl_good(i),i=1,nclstr1)
2062           print*,'singles ',(cl_single(i),i=1,nclstr1)           print*,'singlets',(cl_single(i),i=1,nclstr1)
2063           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
2064        endif        endif
2065                
# Line 1598  c      double precision xm3,ym3,zm3 Line 2119  c      double precision xm3,ym3,zm3
2119        real xc,zc,radius        real xc,zc,radius
2120  *     -----------------------------  *     -----------------------------
2121    
2122          if(DEBUG.EQ.1)print*,'cp_to_doubtrip:'
2123    
2124  *     --------------------------------------------  *     --------------------------------------------
2125  *     put a limit to the maximum number of couples  *     put a limit to the maximum number of couples
# Line 1606  c      double precision xm3,ym3,zm3 Line 2128  c      double precision xm3,ym3,zm3
2128  *     --------------------------------------------  *     --------------------------------------------
2129        do ip=1,nplanes        do ip=1,nplanes
2130           if(ncp_plane(ip).gt.ncouplelimit)then           if(ncp_plane(ip).gt.ncouplelimit)then
2131              mask_view(nviewx(ip)) = 8  c            mask_view(nviewx(ip)) = 8
2132              mask_view(nviewy(ip)) = 8  c            mask_view(nviewy(ip)) = 8
2133                mask_view(nviewx(ip)) = mask_view(nviewx(ip)) + 2**7
2134                mask_view(nviewy(ip)) = mask_view(nviewy(ip)) + 2**7
2135           endif           endif
2136        enddo        enddo
2137    
# Line 1623  c      double precision xm3,ym3,zm3 Line 2147  c      double precision xm3,ym3,zm3
2147                 icx1=clx(ip1,icp1)                 icx1=clx(ip1,icp1)
2148                 icy1=cly(ip1,icp1)                 icy1=cly(ip1,icp1)
2149  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)
2150                 call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)  c               call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)
2151                   call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.)
2152                 xm1=xPAM                 xm1=xPAM
2153                 ym1=yPAM                 ym1=yPAM
2154                 zm1=zPAM                                   zm1=zPAM                  
# Line 1632  c     print*,'***',is1,xm1,ym1,zm1 Line 2157  c     print*,'***',is1,xm1,ym1,zm1
2157                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
2158                    if(  mask_view(nviewx(ip2)).ne.0 .or.                    if(  mask_view(nviewx(ip2)).ne.0 .or.
2159       $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane       $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane
2160                    do is2=1,2    !loop on sensors -ndblt COPPIA 2                    do is2=1,2    !loop on sensors -ndblt COPPIA 2                    
                       
2161                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
2162                          icx2=clx(ip2,icp2)                          icx2=clx(ip2,icp2)
2163                          icy2=cly(ip2,icp2)                          icy2=cly(ip2,icp2)
2164  c                        call xyz_PAM  c                        call xyz_PAM
2165  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)
2166    c                        call xyz_PAM
2167    c     $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)
2168                          call xyz_PAM                          call xyz_PAM
2169       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.)
2170                          xm2=xPAM                          xm2=xPAM
2171                          ym2=yPAM                          ym2=yPAM
2172                          zm2=zPAM                          zm2=zPAM
# Line 1650  c     $                       (icx2,icy2 Line 2176  c     $                       (icx2,icy2
2176  *     (2 couples needed)  *     (2 couples needed)
2177  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2178                          if(ndblt.eq.ndblt_max)then                          if(ndblt.eq.ndblt_max)then
2179                             if(verbose)print*,                             if(verbose.eq.1)print*,
2180       $                          '** warning ** number of identified '//       $                          '** warning ** number of identified '//
2181       $                          'doublets exceeds vector dimention '       $                          'doublets exceeds vector dimention '
2182       $                          ,'( ',ndblt_max,' )'       $                          ,'( ',ndblt_max,' )'
2183  c                           good2=.false.  c                           good2=.false.
2184  c                           goto 880 !fill ntp and go to next event  c                           goto 880 !fill ntp and go to next event
2185                             do iv=1,12                             do iv=1,12
2186                                mask_view(iv) = 3  c                              mask_view(iv) = 3
2187                                  mask_view(iv) = mask_view(iv)+ 2**2
2188                             enddo                             enddo
2189                             iflag=1                             iflag=1
2190                             return                             return
# Line 1704  c$$$ Line 2231  c$$$
2231                                   icy3=cly(ip3,icp3)                                   icy3=cly(ip3,icp3)
2232  c                                 call xyz_PAM  c                                 call xyz_PAM
2233  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)
2234    c                                 call xyz_PAM
2235    c     $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)
2236                                   call xyz_PAM                                   call xyz_PAM
2237       $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)       $                                (icx3,icy3,is3,PFAdef,PFAdef
2238         $                                ,0.,0.,0.,0.)
2239                                   xm3=xPAM                                   xm3=xPAM
2240                                   ym3=yPAM                                   ym3=yPAM
2241                                   zm3=zPAM                                   zm3=zPAM
# Line 1726  c     $                                 Line 2256  c     $                                
2256  *     (3 couples needed)  *     (3 couples needed)
2257  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2258                                   if(ntrpt.eq.ntrpt_max)then                                   if(ntrpt.eq.ntrpt_max)then
2259                                      if(verbose)print*,                                      if(verbose.eq.1)print*,
2260       $                     '** warning ** number of identified '//       $                     '** warning ** number of identified '//
2261       $                     'triplets exceeds vector dimention '       $                     'triplets exceeds vector dimention '
2262       $                    ,'( ',ntrpt_max,' )'       $                    ,'( ',ntrpt_max,' )'
2263  c                                    good2=.false.  c                                    good2=.false.
2264  c                                    goto 880 !fill ntp and go to next event  c                                    goto 880 !fill ntp and go to next event
2265                                      do iv=1,nviews                                      do iv=1,nviews
2266                                         mask_view(iv) = 4  c                                       mask_view(iv) = 4
2267                                           mask_view(iv)=mask_view(iv)+ 2**3
2268                                      enddo                                      enddo
2269                                      iflag=1                                      iflag=1
2270                                      return                                      return
# Line 1787  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp Line 2318  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp
2318   10   continue   10   continue
2319        enddo                     !end loop on planes  - COPPIA 1        enddo                     !end loop on planes  - COPPIA 1
2320                
2321        if(DEBUG)then        if(DEBUG.EQ.1)then
2322           print*,'--- doublets ',ndblt           print*,'--- doublets ',ndblt
2323           print*,'--- triplets ',ntrpt           print*,'--- triplets ',ntrpt
2324        endif        endif
# Line 1834  c      include 'momanhough_init.f' Line 2365  c      include 'momanhough_init.f'
2365        integer cp_useds1(ncouplemaxtot) ! sensor 1        integer cp_useds1(ncouplemaxtot) ! sensor 1
2366        integer cp_useds2(ncouplemaxtot) ! sensor 2        integer cp_useds2(ncouplemaxtot) ! sensor 2
2367    
2368          if(DEBUG.EQ.1)print*,'doub_to_YZcloud:'
2369    
2370  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2371  *     classification of DOUBLETS  *     classification of DOUBLETS
# Line 1960  c         if(ncpused.lt.ncpyz_min)goto 2 Line 2492  c         if(ncpused.lt.ncpyz_min)goto 2
2492  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2493    
2494           if(nclouds_yz.ge.ncloyz_max)then           if(nclouds_yz.ge.ncloyz_max)then
2495              if(verbose)print*,              if(verbose.eq.1)print*,
2496       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2497       $           'YZ clouds exceeds vector dimention '       $           'YZ clouds exceeds vector dimention '
2498       $           ,'( ',ncloyz_max,' )'       $           ,'( ',ncloyz_max,' )'
2499  c               good2=.false.  c               good2=.false.
2500  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2501              do iv=1,nviews              do iv=1,nviews
2502                 mask_view(iv) = 5  c               mask_view(iv) = 5
2503                   mask_view(iv) = mask_view(iv) + 2**4
2504              enddo              enddo
2505              iflag=1              iflag=1
2506              return              return
# Line 1987  c     ptcloud_yz_nt(nclouds_yz)=npt Line 2520  c     ptcloud_yz_nt(nclouds_yz)=npt
2520  c     print*,'>> ',ipt,db_all(ipt)  c     print*,'>> ',ipt,db_all(ipt)
2521           enddo             enddo  
2522           npt_tot=npt_tot+npt           npt_tot=npt_tot+npt
2523           if(DEBUG)then           if(DEBUG.EQ.1)then
2524              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
2525              print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'              print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'
2526              print*,'- alfayz1 ',alfayz1_av(nclouds_yz)              print*,'- alfayz1  ',alfayz1_av(nclouds_yz)
2527              print*,'- alfayz2 ',alfayz2_av(nclouds_yz)              print*,'- alfayz2  ',alfayz2_av(nclouds_yz)
2528              print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)              print*,'cp_useds1  ',(cp_useds1(icp),icp=1,ncp_tot)
2529              print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)              print*,'cp_useds2  ',(cp_useds2(icp),icp=1,ncp_tot)
2530              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)              print*,'cpcloud_yz '
2531         $           ,(cpcloud_yz(nclouds_yz,icp),icp=1,ncp_tot)
2532                print*,'hit_plane  ',(hit_plane(ip),ip=1,nplanes)
2533  c$$$            print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = '  c$$$            print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = '
2534  c$$$     $           ,ptcloud_yz(nclouds_yz)  c$$$     $           ,ptcloud_yz(nclouds_yz)
2535  c$$$            print*,'nt-uple: db_cloud(...) = '  c$$$            print*,'nt-uple: db_cloud(...) = '
# Line 2012  c$$$     $           ,(db_cloud(iii),iii Line 2547  c$$$     $           ,(db_cloud(iii),iii
2547          goto 90                          goto 90                
2548        endif                            endif                    
2549                
2550        if(DEBUG)then        if(DEBUG.EQ.1)then
2551           print*,'---------------------- '           print*,'---------------------- '
2552           print*,'Y-Z total clouds ',nclouds_yz           print*,'Y-Z total clouds ',nclouds_yz
2553           print*,' '           print*,' '
# Line 2061  c      include 'momanhough_init.f' Line 2596  c      include 'momanhough_init.f'
2596        integer cp_useds1(ncouplemaxtot) ! sensor 1        integer cp_useds1(ncouplemaxtot) ! sensor 1
2597        integer cp_useds2(ncouplemaxtot) ! sensor 2        integer cp_useds2(ncouplemaxtot) ! sensor 2
2598    
2599          if(DEBUG.EQ.1)print*,'trip_to_XZcloud:'
2600    
2601  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2602  *     classification of TRIPLETS  *     classification of TRIPLETS
2603  *     according to distance in parameter space  *     according to distance in parameter space
# Line 2118  c         tr_temp(1)=itr1 Line 2655  c         tr_temp(1)=itr1
2655       $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2                     $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2              
2656                 distance = sqrt(distance)                 distance = sqrt(distance)
2657                                
2658                 if(distance.lt.cutdistxz)then  *     ------------------------------------------------------------------------
2659    *     FORCE INCLUSION OF TRIPLETS COMPOSED BY SAME COUPLES, IGNORING THE IMAGE
2660    *     ------------------------------------------------------------------------
2661    *     (added in august 2007)
2662                   istrimage=0
2663                   if(
2664         $              abs(cpxz1(itrref)).eq.abs(cpxz1(itr2)).and.
2665         $              abs(cpxz2(itrref)).eq.abs(cpxz2(itr2)).and.
2666         $              abs(cpxz3(itrref)).eq.abs(cpxz3(itr2)).and.
2667         $              .true.)istrimage=1
2668    
2669                   if(distance.lt.cutdistxz.or.istrimage.eq.1)then
2670  c     print*,idb1,idb2,distance,' cloud ',nclouds_yz  c     print*,idb1,idb2,distance,' cloud ',nclouds_yz
2671                    if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1                    if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1
2672                    if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1                    if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1
# Line 2177  c     print*,'check cp_used' Line 2725  c     print*,'check cp_used'
2725           enddo           enddo
2726  c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet  c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet
2727           if(npt.lt.nptxz_min)goto 22288     !next triplet           if(npt.lt.nptxz_min)goto 22288     !next triplet
2728           if(nplused.lt.nplxz_min)goto 22288 !next doublet           if(nplused.lt.nplxz_min)goto 22288 !next triplet
2729                    
2730  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2731  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2732           if(nclouds_xz.ge.ncloxz_max)then           if(nclouds_xz.ge.ncloxz_max)then
2733              if(verbose)print*,              if(verbose.eq.1)print*,
2734       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2735       $           'XZ clouds exceeds vector dimention '       $           'XZ clouds exceeds vector dimention '
2736       $           ,'( ',ncloxz_max,' )'       $           ,'( ',ncloxz_max,' )'
2737  c     good2=.false.  c     good2=.false.
2738  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2739              do iv=1,nviews              do iv=1,nviews
2740                 mask_view(iv) = 6  c               mask_view(iv) = 6
2741                   mask_view(iv) =  mask_view(iv) + 2**5
2742              enddo              enddo
2743              iflag=1              iflag=1
2744              return              return
# Line 2208  c     goto 880         !fill ntp and go Line 2757  c     goto 880         !fill ntp and go
2757           enddo           enddo
2758           npt_tot=npt_tot+npt           npt_tot=npt_tot+npt
2759                    
2760           if(DEBUG)then           if(DEBUG.EQ.1)then
2761              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
2762              print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'                            print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'              
2763              print*,'- alfaxz1 ',alfaxz1_av(nclouds_xz)              print*,'- alfaxz1  ',alfaxz1_av(nclouds_xz)
2764              print*,'- alfaxz2 ',alfaxz2_av(nclouds_xz)              print*,'- alfaxz2  ',alfaxz2_av(nclouds_xz)
2765              print*,'- alfaxz3 ',alfaxz3_av(nclouds_xz)              print*,'- alfaxz3  ',alfaxz3_av(nclouds_xz)
2766              print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)              print*,'cp_useds1  ',(cp_useds1(icp),icp=1,ncp_tot)
2767              print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)              print*,'cp_useds2  ',(cp_useds2(icp),icp=1,ncp_tot)
2768                print*,'cpcloud_xz '
2769         $           ,(cpcloud_xz(nclouds_xz,icp),icp=1,ncp_tot)
2770              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)
2771  c$$$            print*,'nt-uple: ptcloud_xz(',nclouds_xz,') = '  c$$$            print*,'nt-uple: ptcloud_xz(',nclouds_xz,') = '
2772  c$$$     $           ,ptcloud_xz(nclouds_xz)  c$$$     $           ,ptcloud_xz(nclouds_xz)
# Line 2233  c$$$     $           ,(tr_cloud(iii),iii Line 2784  c$$$     $           ,(tr_cloud(iii),iii
2784           goto 91                           goto 91                
2785         endif                             endif                    
2786                
2787        if(DEBUG)then        if(DEBUG.EQ.1)then
2788           print*,'---------------------- '           print*,'---------------------- '
2789           print*,'X-Z total clouds ',nclouds_xz           print*,'X-Z total clouds ',nclouds_xz
2790           print*,' '           print*,' '
# Line 2254  c$$$     $           ,(tr_cloud(iii),iii Line 2805  c$$$     $           ,(tr_cloud(iii),iii
2805  **************************************************  **************************************************
2806    
2807        subroutine clouds_to_ctrack(iflag)        subroutine clouds_to_ctrack(iflag)
 c*****************************************************  
 c     02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
2808    
2809        include 'commontracker.f'        include 'commontracker.f'
2810        include 'level1.f'        include 'level1.f'
# Line 2264  c*************************************** Line 2812  c***************************************
2812        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2813        include 'common_mini_2.f'        include 'common_mini_2.f'
2814        include 'common_mech.f'        include 'common_mech.f'
2815  c      include 'momanhough_init.f'  
2816    
2817    
2818  *     output flag  *     output flag
# Line 2288  c      include 'momanhough_init.f' Line 2836  c      include 'momanhough_init.f'
2836  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2837  *     variables for track fitting  *     variables for track fitting
2838        double precision AL_INI(5)        double precision AL_INI(5)
 c      double precision tath  
2839  *     -----------------------------------------------------------  *     -----------------------------------------------------------
 c      real fitz(nplanes)        !z coordinates of the planes in cm  
2840    
2841          if(DEBUG.EQ.1)print*,'clouds_to_ctrack:'
2842    
2843    
2844        ntracks=0                 !counter of track candidates        ntracks=0                 !counter of track candidates
# Line 2313  c      real fitz(nplanes)        !z coor Line 2860  c      real fitz(nplanes)        !z coor
2860                 enddo                 enddo
2861              enddo              enddo
2862              ncp_ok=0              ncp_ok=0
2863              do icp=1,ncp_tot    !loop on couples              do icp=1,ncp_tot    !loop over couples
2864  *     get info on  *     get info on
2865                 cpintersec(icp)=min(                 cpintersec(icp)=min(
2866       $              cpcloud_yz(iyz,icp),       $              cpcloud_yz(iyz,icp),
# Line 2322  c      real fitz(nplanes)        !z coor Line 2869  c      real fitz(nplanes)        !z coor
2869       $    (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.
2870       $    (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.
2871       $              .false.)cpintersec(icp)=0       $              .false.)cpintersec(icp)=0
2872    *     cpintersec is >0 if yz and xz clouds contain the same image of couple icp
2873                 if(cpintersec(icp).ne.0)then                 if(cpintersec(icp).ne.0)then
2874                    ncp_ok=ncp_ok+1                      ncp_ok=ncp_ok+1  
2875                                        
# Line 2354  c      real fitz(nplanes)        !z coor Line 2902  c      real fitz(nplanes)        !z coor
2902                 nplused=nplused+ hit_plane(ip)                 nplused=nplused+ hit_plane(ip)
2903              enddo              enddo
2904                            
 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  
2905                            
2906              if(DEBUG)then              if(DEBUG.EQ.1)then
2907                 print*,'Combination ',iyz,ixz                 print*,'Combination ',iyz,ixz
2908       $              ,' db ',ptcloud_yz(iyz)       $              ,' db ',ptcloud_yz(iyz)
2909       $              ,' tr ',ptcloud_xz(ixz)       $              ,' tr ',ptcloud_xz(ixz)
2910       $              ,'  -----> # matching couples ',ncp_ok       $              ,'  -----> # matching couples ',ncp_ok
2911              endif              endif
2912    
2913    c            if(nplused.lt.nplxz_min)goto 888 !next combination
2914                if(nplused.lt.nplyz_min)goto 888 !next combination
2915                if(ncp_ok.lt.ncpok)goto 888 !next combination
2916    
2917  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
2918  c$$$  print*,'Configurazione cluster XZ'  c$$$  print*,'Configurazione cluster XZ'
2919  c$$$  print*,'1 -- ',(clx(1,i),i=1,ncp_plane(1))  c$$$  print*,'1 -- ',(clx(1,i),i=1,ncp_plane(1))
# Line 2409  c$$$            AL_INI(5) = (1.e2*alfaxz Line 2959  c$$$            AL_INI(5) = (1.e2*alfaxz
2959  c$$$              c$$$            
2960  c$$$            if(AL_INI(5).gt.defmax)goto 888 !next cloud  c$$$            if(AL_INI(5).gt.defmax)goto 888 !next cloud
2961                                                    
2962              if(DEBUG)then              if(DEBUG.EQ.1)then
2963                   print*,'track candidate', ntracks+1
2964                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))
2965                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))
2966                 print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4))                 print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4))
# Line 2442  c$$$            if(AL_INI(5).gt.defmax)g Line 2993  c$$$            if(AL_INI(5).gt.defmax)g
2993                                hit_plane(6)=icp6                                hit_plane(6)=icp6
2994                                if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6                                if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6
2995                                                                
2996                                  *                             ---------------------------------------
2997    *                             check if this group of couples has been
2998    *                             already fitted    
2999    *                             ---------------------------------------
3000                                  do ica=1,ntracks
3001                                     isthesame=1
3002                                     do ip=1,NPLANES
3003                                        if(hit_plane(ip).ne.0)then
3004                                           if(  CP_STORE(nplanes-ip+1,ica)
3005         $                                      .ne.
3006         $                                      cp_match(ip,hit_plane(ip)) )
3007         $                                      isthesame=0
3008                                        else
3009                                           if(  CP_STORE(nplanes-ip+1,ica)
3010         $                                      .ne.
3011         $                                      0 )
3012         $                                      isthesame=0
3013                                        endif
3014                                     enddo
3015                                     if(isthesame.eq.1)then
3016                                        if(DEBUG.eq.1)
3017         $                                   print*,'(already fitted)'
3018                                        goto 666 !jump to next combination
3019                                     endif
3020                                  enddo
3021    
3022                                call track_init !init TRACK common                                call track_init !init TRACK common
3023    
3024                                do ip=1,nplanes !loop on planes                                do ip=1,nplanes !loop on planes (bottom to top)
3025                                   if(hit_plane(ip).ne.0)then                                   if(hit_plane(ip).ne.0)then
3026                                      id=cp_match(ip,hit_plane(ip))                                      id=cp_match(ip,hit_plane(ip))
3027                                      is=is_cp(id)                                      is=is_cp(id)
# Line 2459  c$$$            if(AL_INI(5).gt.defmax)g Line 3035  c$$$            if(AL_INI(5).gt.defmax)g
3035  *                                   *************************  *                                   *************************
3036  c                                    call xyz_PAM(icx,icy,is,  c                                    call xyz_PAM(icx,icy,is,
3037  c     $                                   'COG2','COG2',0.,0.)  c     $                                   'COG2','COG2',0.,0.)
3038    c                                    call xyz_PAM(icx,icy,is, !(1)
3039    c     $                                   PFAdef,PFAdef,0.,0.) !(1)
3040                                      call xyz_PAM(icx,icy,is, !(1)                                      call xyz_PAM(icx,icy,is, !(1)
3041       $                                   PFAdef,PFAdef,0.,0.) !(1)       $                                   PFAdef,PFAdef,0.,0.,0.,0.)
3042  *                                   *************************  *                                   *************************
3043  *                                   -----------------------------  *                                   -----------------------------
3044                                      xgood(nplanes-ip+1)=1.                                      xgood(nplanes-ip+1)=1.
# Line 2487  c$$$                              enddo Line 3065  c$$$                              enddo
3065                                ifail=0 !error flag in chi^2 computation                                ifail=0 !error flag in chi^2 computation
3066                                jstep=0 !number of  minimization steps                                jstep=0 !number of  minimization steps
3067                                iprint=0                                iprint=0
3068  c                              if(DEBUG)iprint=1  c                              if(DEBUG.EQ.1)iprint=1
3069                                if(DEBUG)iprint=2                                if(DEBUG.EQ.1)iprint=2
3070                                call mini2(jstep,ifail,iprint)                                call mini2(jstep,ifail,iprint)
3071                                if(ifail.ne.0) then                                if(ifail.ne.0) then
3072                                   if(DEBUG)then                                   if(DEBUG.EQ.1)then
3073                                      print *,                                      print *,
3074       $                              '*** MINIMIZATION FAILURE *** '       $                              '*** MINIMIZATION FAILURE *** '
3075       $                              //'(clouds_to_ctrack)'       $                              //'(clouds_to_ctrack)'
# Line 2516  c                                 chi2=- Line 3094  c                                 chi2=-
3094  *     --------------------------  *     --------------------------
3095                                if(ntracks.eq.NTRACKSMAX)then                                if(ntracks.eq.NTRACKSMAX)then
3096                                                                    
3097                                   if(verbose)print*,                                   if(verbose.eq.1)print*,
3098       $                 '** warning ** number of candidate tracks '//       $                 '** warning ** number of candidate tracks '//
3099       $                 ' exceeds vector dimension '       $                 ' exceeds vector dimension '
3100       $                ,'( ',NTRACKSMAX,' )'       $                ,'( ',NTRACKSMAX,' )'
3101  c                                 good2=.false.  c                                 good2=.false.
3102  c                                 goto 880 !fill ntp and go to next event                      c                                 goto 880 !fill ntp and go to next event                    
3103                                   do iv=1,nviews                                   do iv=1,nviews
3104                                      mask_view(iv) = 7  c                                    mask_view(iv) = 7
3105                                        mask_view(iv) = mask_view(iv) + 2**6
3106                                   enddo                                   enddo
3107                                   iflag=1                                   iflag=1
3108                                   return                                   return
# Line 2531  c                                 goto 8 Line 3110  c                                 goto 8
3110                                                                
3111                                ntracks = ntracks + 1                                ntracks = ntracks + 1
3112                                                                
3113  c$$$                              ndof=0                                                                do ip=1,nplanes !top to bottom
3114                                do ip=1,nplanes  
 c$$$                                 ndof=ndof  
 c$$$     $                                +int(xgood(ip))  
 c$$$     $                                +int(ygood(ip))  
3115                                   XV_STORE(ip,ntracks)=sngl(xv(ip))                                   XV_STORE(ip,ntracks)=sngl(xv(ip))
3116                                   YV_STORE(ip,ntracks)=sngl(yv(ip))                                   YV_STORE(ip,ntracks)=sngl(yv(ip))
3117                                   ZV_STORE(ip,ntracks)=sngl(zv(ip))                                                                       ZV_STORE(ip,ntracks)=sngl(zv(ip))                                    
# Line 2551  c$$$     $                               Line 3127  c$$$     $                              
3127                                   AYV_STORE(ip,ntracks)=sngl(ayv(ip))                                   AYV_STORE(ip,ntracks)=sngl(ayv(ip))
3128                                   XGOOD_STORE(ip,ntracks)=sngl(xgood(ip))                                   XGOOD_STORE(ip,ntracks)=sngl(xgood(ip))
3129                                   YGOOD_STORE(ip,ntracks)=sngl(ygood(ip))                                   YGOOD_STORE(ip,ntracks)=sngl(ygood(ip))
3130    *                                NB! hit_plane is defined from bottom to top
3131                                   if(hit_plane(ip).ne.0)then                                   if(hit_plane(ip).ne.0)then
3132                                      CP_STORE(nplanes-ip+1,ntracks)=                                      CP_STORE(nplanes-ip+1,ntracks)=
3133       $                                   cp_match(ip,hit_plane(ip))       $                                   cp_match(ip,hit_plane(ip))
3134                                        SENSOR_STORE(nplanes-ip+1,ntracks)
3135         $                              = is_cp(cp_match(ip,hit_plane(ip)))
3136                                        LADDER_STORE(nplanes-ip+1,ntracks)
3137         $                                   = LADDER(
3138         $                                   clx(ip,icp_cp(
3139         $                                   cp_match(ip,hit_plane(ip)
3140         $                                   ))));
3141                                   else                                   else
3142                                      CP_STORE(nplanes-ip+1,ntracks)=0                                      CP_STORE(nplanes-ip+1,ntracks)=0
3143                                        SENSOR_STORE(nplanes-ip+1,ntracks)=0
3144                                        LADDER_STORE(nplanes-ip+1,ntracks)=0
3145                                   endif                                   endif
3146                                   CLS_STORE(nplanes-ip+1,ntracks)=0                                   BX_STORE(ip,ntracks)=0!I dont need it now
3147                                     BY_STORE(ip,ntracks)=0!I dont need it now
3148                                     CLS_STORE(ip,ntracks)=0
3149                                   do i=1,5                                   do i=1,5
3150                                      AL_STORE(i,ntracks)=sngl(AL(i))                                      AL_STORE(i,ntracks)=sngl(AL(i))
3151                                   enddo                                   enddo
3152                                enddo                                enddo
3153                                                                
 c$$$  *                             Number of Degree Of Freedom  
 c$$$  ndof=ndof-5                            
 c$$$  *                             reduced chi^2  
 c$$$  rchi2=chi2/dble(ndof)  
3154                                RCHI2_STORE(ntracks)=chi2                                RCHI2_STORE(ntracks)=chi2
3155                                                                
3156  *     --------------------------------  *     --------------------------------
# Line 2590  c$$$  rchi2=chi2/dble(ndof) Line 3174  c$$$  rchi2=chi2/dble(ndof)
3174           return           return
3175        endif        endif
3176                
3177        if(DEBUG)then  c$$$      if(DEBUG.EQ.1)then
3178           print*,'****** TRACK CANDIDATES ***********'  c$$$         print*,'****** TRACK CANDIDATES ***********'
3179           print*,'#         R. chi2        RIG'  c$$$         print*,'#         R. chi2        RIG'
3180           do i=1,ntracks  c$$$         do i=1,ntracks
3181              print*,i,' --- ',rchi2_store(i),' --- '  c$$$            print*,i,' --- ',rchi2_store(i),' --- '
3182       $           ,1./abs(AL_STORE(5,i))  c$$$     $           ,1./abs(AL_STORE(5,i))
3183           enddo  c$$$         enddo
3184           print*,'***********************************'  c$$$         print*,'***********************************'
3185    c$$$      endif
3186          if(DEBUG.EQ.1)then
3187            print*,'****** TRACK CANDIDATES *****************'
3188            print*,'#         R. chi2        RIG         ndof'
3189            do i=1,ntracks
3190              ndof=0                !(1)
3191              do ii=1,nplanes       !(1)
3192                ndof=ndof           !(1)
3193         $           +int(xgood_store(ii,i)) !(1)
3194         $           +int(ygood_store(ii,i)) !(1)
3195              enddo                 !(1)
3196              print*,i,' --- ',rchi2_store(i),' --- '
3197         $         ,1./abs(AL_STORE(5,i)),' --- ',ndof
3198            enddo
3199            print*,'*****************************************'
3200        endif        endif
3201                
3202                
# Line 2616  c$$$  rchi2=chi2/dble(ndof) Line 3215  c$$$  rchi2=chi2/dble(ndof)
3215    
3216        subroutine refine_track(ibest)        subroutine refine_track(ibest)
3217    
 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******************************************************  
3218    
3219        include 'commontracker.f'        include 'commontracker.f'
3220        include 'level1.f'        include 'level1.f'
# Line 2628  c*************************************** Line 3222  c***************************************
3222        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3223        include 'common_mini_2.f'        include 'common_mini_2.f'
3224        include 'common_mech.f'        include 'common_mech.f'
 c      include 'momanhough_init.f'  
 c      include 'level1.f'  
3225        include 'calib.f'        include 'calib.f'
3226    
   
3227  *     flag to chose PFA  *     flag to chose PFA
3228        character*10 PFA        character*10 PFA
3229        common/FINALPFA/PFA        common/FINALPFA/PFA
3230    
3231          real k(6)
3232          DATA k/1.099730,0.418900,0.220939,0.220907,0.418771,1.100674/
3233    
3234          real xp,yp,zp
3235          real xyzp(3),bxyz(3)
3236          equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3))
3237    
3238          if(DEBUG.EQ.1)print*,'refine_track:'
3239  *     =================================================  *     =================================================
3240  *     new estimate of positions using ETA algorithm  *     new estimate of positions using ETA algorithm
3241  *                          and  *                          and
# Line 2645  c      include 'level1.f' Line 3244  c      include 'level1.f'
3244        call track_init        call track_init
3245        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3246    
3247             xP=XV_STORE(nplanes-ip+1,ibest)
3248             yP=YV_STORE(nplanes-ip+1,ibest)
3249             zP=ZV_STORE(nplanes-ip+1,ibest)
3250             call gufld(xyzp,bxyz)
3251             BX_STORE(nplanes-ip+1,ibest)=bxyz(1)
3252             BY_STORE(nplanes-ip+1,ibest)=bxyz(2)
3253    c$$$  bxyz(1)=0
3254    c$$$         bxyz(2)=0
3255    c$$$         bxyz(3)=0
3256  *     |||||||||||||||||||||||||||||||||||||||||||||||||  *     |||||||||||||||||||||||||||||||||||||||||||||||||
3257  *     -------------------------------------------------  *     -------------------------------------------------
3258  *     If the plane has been already included, it just  *     If the plane has been already included, it just
# Line 2665  c      include 'level1.f' Line 3273  c      include 'level1.f'
3273       $           ,ip_cp(id),ip       $           ,ip_cp(id),ip
3274              icx=clx(ip,icp)              icx=clx(ip,icp)
3275              icy=cly(ip,icp)              icy=cly(ip,icp)
3276    c            call xyz_PAM(icx,icy,is,
3277    c     $           PFA,PFA,
3278    c     $           AXV_STORE(nplanes-ip+1,ibest),
3279    c     $           AYV_STORE(nplanes-ip+1,ibest))
3280              call xyz_PAM(icx,icy,is,              call xyz_PAM(icx,icy,is,
 c     $           'ETA2','ETA2',  
3281       $           PFA,PFA,       $           PFA,PFA,
3282       $           AXV_STORE(nplanes-ip+1,ibest),       $           AXV_STORE(nplanes-ip+1,ibest),
3283       $           AYV_STORE(nplanes-ip+1,ibest))       $           AYV_STORE(nplanes-ip+1,ibest),
3284  c$$$  call xyz_PAM(icx,icy,is,       $           bxyz(1),
3285  c$$$  $              'COG2','COG2',       $           bxyz(2)
3286  c$$$  $              0.,       $           )
3287  c$$$  $              0.)  
3288              xm(nplanes-ip+1) = xPAM              xm(nplanes-ip+1) = xPAM
3289              ym(nplanes-ip+1) = yPAM              ym(nplanes-ip+1) = yPAM
3290              zm(nplanes-ip+1) = zPAM              zm(nplanes-ip+1) = zPAM
# Line 2682  c$$$  $              0.) Line 3293  c$$$  $              0.)
3293              resx(nplanes-ip+1) = resxPAM              resx(nplanes-ip+1) = resxPAM
3294              resy(nplanes-ip+1) = resyPAM              resy(nplanes-ip+1) = resyPAM
3295    
3296  c            dedxtrk(nplanes-ip+1) = (dedx(icx)+dedx(icy))/2. !(1)              dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx))
3297              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)  
3298                            
3299  *     |||||||||||||||||||||||||||||||||||||||||||||||||  *     |||||||||||||||||||||||||||||||||||||||||||||||||
3300  *     -------------------------------------------------  *     -------------------------------------------------
# Line 2699  c            dedxtrk(nplanes-ip+1) = (de Line 3309  c            dedxtrk(nplanes-ip+1) = (de
3309                                
3310  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3311  *     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)  
3312              call whichsensor(ip,xP,yP,nldt,ist)              call whichsensor(ip,xP,yP,nldt,ist)
3313  *     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
3314              if(nldt.eq.0.or.ist.eq.0)goto 133              if(nldt.eq.0.or.ist.eq.0)goto 133
3315    
3316                SENSOR_STORE(nplanes-ip+1,IBEST)=ist
3317                LADDER_STORE(nplanes-ip+1,IBEST)=nldt
3318  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3319    
3320              if(DEBUG)then              if(DEBUG.EQ.1)then
3321                 print*,                 print*,
3322       $              '------ Plane ',ip,' intersected on LADDER ',nldt       $              '------ Plane ',ip,' intersected on LADDER ',nldt
3323       $              ,' SENSOR ',ist       $              ,' SENSOR ',ist
# Line 2718  c            dedxtrk(nplanes-ip+1) = (de Line 3328  c            dedxtrk(nplanes-ip+1) = (de
3328  *     ===========================================  *     ===========================================
3329  *     STEP 1 >>>>>>>  try to include a new couple  *     STEP 1 >>>>>>>  try to include a new couple
3330  *     ===========================================  *     ===========================================
3331  c            if(DEBUG)print*,'>>>> try to include a new couple'  c            if(DEBUG.EQ.1)print*,'>>>> try to include a new couple'
3332              distmin=1000000.              distmin=1000000.
3333              xmm = 0.              xmm = 0.
3334              ymm = 0.              ymm = 0.
# Line 2741  c     $              cl_used(icy).eq.1.o Line 3351  c     $              cl_used(icy).eq.1.o
3351  *            *          
3352                 call xyz_PAM(icx,icy,ist,                 call xyz_PAM(icx,icy,ist,
3353       $              PFA,PFA,       $              PFA,PFA,
 c     $              'ETA2','ETA2',  
3354       $              AXV_STORE(nplanes-ip+1,ibest),       $              AXV_STORE(nplanes-ip+1,ibest),
3355       $              AYV_STORE(nplanes-ip+1,ibest))       $              AYV_STORE(nplanes-ip+1,ibest),
3356         $              bxyz(1),
3357         $              bxyz(2)
3358         $              )
3359                                
3360                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3361                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3362                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
3363                 if(DEBUG)print*,'( couple ',id                 if(DEBUG.EQ.1)print*,'( couple ',id
3364       $              ,' ) normalized distance ',distance       $              ,' ) distance ',distance
3365                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3366                    xmm = xPAM                    xmm = xPAM
3367                    ymm = yPAM                    ymm = yPAM
# Line 2758  c     $              'ETA2','ETA2', Line 3370  c     $              'ETA2','ETA2',
3370                    rymm = resyPAM                    rymm = resyPAM
3371                    distmin = distance                    distmin = distance
3372                    idm = id                                      idm = id                  
3373  c                 dedxmm = (dedx(icx)+dedx(icy))/2. !(1)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3374                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3375                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)  c     QUIQUI --> non devo moltiplicare per clinc?!?!?!
3376                      clincnewc=10*sqrt(rymm**2+rxmm**2 !QUIQUI
3377         $                 +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2))) !QUIQUI
3378                 endif                 endif
3379   1188          continue   1188          continue
3380              enddo               !end loop on couples on plane icp              enddo               !end loop on couples on plane icp
3381              if(distmin.le.clinc)then                    c            if(distmin.le.clinc)then     !QUIQUI              
3382                if(distmin.le.clincnewc)then     !QUIQUI              
3383  *              -----------------------------------  *              -----------------------------------
3384                 xm(nplanes-ip+1) = xmm         !<<<                 xm(nplanes-ip+1) = xmm !<<<
3385                 ym(nplanes-ip+1) = ymm         !<<<                 ym(nplanes-ip+1) = ymm !<<<
3386                 zm(nplanes-ip+1) = zmm         !<<<                 zm(nplanes-ip+1) = zmm !<<<
3387                 xgood(nplanes-ip+1) = 1        !<<<                 xgood(nplanes-ip+1) = 1 !<<<
3388                 ygood(nplanes-ip+1) = 1        !<<<                 ygood(nplanes-ip+1) = 1 !<<<
3389                 resx(nplanes-ip+1)=rxmm        !<<<                 resx(nplanes-ip+1)=rxmm !<<<
3390                 resy(nplanes-ip+1)=rymm        !<<<                 resy(nplanes-ip+1)=rymm !<<<
3391  c              dedxtrk(nplanes-ip+1) = dedxmm !<<<  !(1)                 dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3392                 dedxtrk_x(nplanes-ip+1) = dedxmmx    !(1)                 dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
                dedxtrk_y(nplanes-ip+1) = dedxmmy    !(1)  
3393  *              -----------------------------------  *              -----------------------------------
3394                 CP_STORE(nplanes-ip+1,ibest)=idm                       CP_STORE(nplanes-ip+1,ibest)=idm      
3395                 if(DEBUG)print*,'%%%% included couple ',idm                 if(DEBUG.EQ.1)print*,'%%%% included couple ',idm
3396       $              ,' (norm.dist.= ',distmin,', cut ',clinc,' )'       $              ,' (dist.= ',distmin,', cut ',clinc,' )'
3397                 goto 133         !next plane                 goto 133         !next plane
3398              endif              endif
3399  *     ================================================  *     ================================================
3400  *     STEP 2 >>>>>>>  try to include a single cluster  *     STEP 2 >>>>>>>  try to include a single cluster
3401  *                     either from a couple or single  *                     either from a couple or single
3402  *     ================================================  *     ================================================
3403  c            if(DEBUG)print*,'>>>> try to include a new cluster'  c            if(DEBUG.EQ.1)print*,'>>>> try to include a new cluster'
3404              distmin=1000000.              distmin=1000000.
3405              xmm_A = 0.          !---------------------------              xmm_A = 0.          !---------------------------
3406              ymm_A = 0.          ! init variables that              ymm_A = 0.          ! init variables that
# Line 2812  c            if(DEBUG)print*,'>>>> try t Line 3426  c            if(DEBUG)print*,'>>>> try t
3426  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
3427                 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)
3428  *                                              !jump to the Y cluster  *                                              !jump to the Y cluster
3429    c               call xyz_PAM(icx,0,ist,
3430    c     $              PFA,PFA,
3431    c     $              AXV_STORE(nplanes-ip+1,ibest),0.)              
3432                 call xyz_PAM(icx,0,ist,                 call xyz_PAM(icx,0,ist,
 c     $              'ETA2','ETA2',  
3433       $              PFA,PFA,       $              PFA,PFA,
3434       $              AXV_STORE(nplanes-ip+1,ibest),0.)                     $              AXV_STORE(nplanes-ip+1,ibest),0.,
3435         $              bxyz(1),
3436         $              bxyz(2)
3437         $              )              
3438                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3439                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3440                 if(DEBUG)print*,'( cl-X ',icx                 if(DEBUG.EQ.1)print*,'( cl-X ',icx
3441       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3442                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3443                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3444                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 2831  c     $              'ETA2','ETA2', Line 3450  c     $              'ETA2','ETA2',
3450                    rymm = resyPAM                    rymm = resyPAM
3451                    distmin = distance                    distmin = distance
3452                    iclm = icx                    iclm = icx
3453  c                  dedxmm = dedx(icx) !(1)  c                  dedxmm = sgnl(icx) !(1)
3454                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3455                    dedxmmy = 0.        !(1)                    dedxmmy = 0.        !(1)
3456                 endif                                   endif                  
3457  11881          continue  11881          continue
# Line 2840  c                  dedxmm = dedx(icx) !( Line 3459  c                  dedxmm = dedx(icx) !(
3459  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
3460                 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)
3461  *                                              !jump to the next couple  *                                              !jump to the next couple
3462    c               call xyz_PAM(0,icy,ist,
3463    c     $              PFA,PFA,
3464    c     $              0.,AYV_STORE(nplanes-ip+1,ibest))
3465                 call xyz_PAM(0,icy,ist,                 call xyz_PAM(0,icy,ist,
 c     $              'ETA2','ETA2',  
3466       $              PFA,PFA,       $              PFA,PFA,
3467       $              0.,AYV_STORE(nplanes-ip+1,ibest))       $              0.,AYV_STORE(nplanes-ip+1,ibest),
3468         $              bxyz(1),
3469         $              bxyz(2)
3470         $              )
3471                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3472                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3473                 if(DEBUG)print*,'( cl-Y ',icy                 if(DEBUG.EQ.1)print*,'( cl-Y ',icy
3474       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3475                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3476                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3477                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 2859  c     $              'ETA2','ETA2', Line 3483  c     $              'ETA2','ETA2',
3483                    rymm = resyPAM                    rymm = resyPAM
3484                    distmin = distance                    distmin = distance
3485                    iclm = icy                    iclm = icy
3486  c                 dedxmm = dedx(icy)  !(1)  c                 dedxmm = sgnl(icy)  !(1)
3487                    dedxmmx = 0.        !(1)                    dedxmmx = 0.        !(1)
3488                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3489                 endif                                   endif                  
3490  11882          continue  11882          continue
3491              enddo               !end loop on cluster inside couples              enddo               !end loop on cluster inside couples
3492  *----- single clusters -----------------------------------------------      *----- single clusters -----------------------------------------------  
3493    c            print*,'## ncls(',ip,') ',ncls(ip)
3494              do ic=1,ncls(ip)    !loop on single clusters              do ic=1,ncls(ip)    !loop on single clusters
3495                 icl=cls(ip,ic)                 icl=cls(ip,ic)
3496  c               if(cl_used(icl).eq.1.or.     !if the cluster is already used  c               if(cl_used(icl).eq.1.or.     !if the cluster is already used
# Line 2874  c               if(cl_used(icl).eq.1.or. Line 3499  c               if(cl_used(icl).eq.1.or.
3499       $              .false.)goto 18882      !jump to the next singlet       $              .false.)goto 18882      !jump to the next singlet
3500                 if(mod(VIEW(icl),2).eq.0)then!<---- X view                 if(mod(VIEW(icl),2).eq.0)then!<---- X view
3501                    call xyz_PAM(icl,0,ist,                    call xyz_PAM(icl,0,ist,
 c     $                 'ETA2','ETA2',  
3502       $                 PFA,PFA,       $                 PFA,PFA,
3503       $                 AXV_STORE(nplanes-ip+1,ibest),0.)       $                 AXV_STORE(nplanes-ip+1,ibest),0.,
3504         $                 bxyz(1),
3505         $                 bxyz(2)
3506         $                 )
3507                 else                         !<---- Y view                 else                         !<---- Y view
3508                    call xyz_PAM(0,icl,ist,                    call xyz_PAM(0,icl,ist,
 c     $                 'ETA2','ETA2',  
3509       $                 PFA,PFA,       $                 PFA,PFA,
3510       $                 0.,AYV_STORE(nplanes-ip+1,ibest))       $                 0.,AYV_STORE(nplanes-ip+1,ibest),
3511         $                 bxyz(1),
3512         $                 bxyz(2)
3513         $                 )
3514                 endif                 endif
3515    
3516                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3517                 distance = distance / RCHI2_STORE(ibest)!<<< MS  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3518                 if(DEBUG)print*,'( cl-s ',icl                 if(DEBUG.EQ.1)print*,'( cl-s ',icl
3519       $              ,' ) normalized distance ',distance       $              ,' ) distance ',distance
3520                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3521    c                  if(DEBUG.EQ.1)print*,'YES'
3522                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3523                    ymm_A = yPAM_A                    ymm_A = yPAM_A
3524                    zmm_A = zPAM_A                    zmm_A = zPAM_A
# Line 2899  c     $                 'ETA2','ETA2', Line 3529  c     $                 'ETA2','ETA2',
3529                    rymm = resyPAM                    rymm = resyPAM
3530                    distmin = distance                      distmin = distance  
3531                    iclm = icl                    iclm = icl
 c                  dedxmm = dedx(icl)                   !(1)  
3532                    if(mod(VIEW(icl),2).eq.0)then !<---- X view                    if(mod(VIEW(icl),2).eq.0)then !<---- X view
3533                       dedxmmx = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmx = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3534                       dedxmmy = 0.                       !(1)                       dedxmmy = 0.                  
3535                    else          !<---- Y view                    else          !<---- Y view
3536                       dedxmmx = 0.                       !(1)                       dedxmmx = 0.                  
3537                       dedxmmy = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmy = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3538                    endif                    endif
3539                 endif                                   endif                  
3540  18882          continue  18882          continue
3541              enddo               !end loop on single clusters              enddo               !end loop on single clusters
3542    c            print*,'## distmin ', distmin,' clinc ',clinc
3543    
3544              if(distmin.le.clinc)then                    c     QUIQUI------------
3545                  c     anche qui: non ci vuole clinc???
3546                 CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<                  if(iclm.ne.0)then
 *              ----------------------------  
 c               print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'  
3547                 if(mod(VIEW(iclm),2).eq.0)then                 if(mod(VIEW(iclm),2).eq.0)then
3548                    XGOOD(nplanes-ip+1)=1.                    clincnew=
3549                    resx(nplanes-ip+1)=rxmm       $                 20*
3550                    if(DEBUG)print*,'%%%% included X-cl ',iclm       $                 sqrt(rxmm**2+RCHI2_STORE(ibest)*k(ip)*cov(1,1))
3551  c                  if(.true.)print*,'%%%% included X-cl ',iclm                 else if(mod(VIEW(iclm),2).ne.0)then
3552       $                 ,'( chi^2, ',RCHI2_STORE(ibest)                    clincnew=
3553       $                 ,', norm.dist.= ',distmin       $                 10*
3554       $                 ,', cut ',clinc,' )'       $                 sqrt(rymm**2+RCHI2_STORE(ibest)*k(ip)*cov(2,2))
3555                 else                 endif
3556                    YGOOD(nplanes-ip+1)=1.  c     QUIQUI------------
3557                    resy(nplanes-ip+1)=rymm                
3558                    if(DEBUG)print*,'%%%% included Y-cl ',iclm                 if(distmin.le.clincnew)then   !QUIQUI
3559  c                  if(.true.)print*,'%%%% included Y-cl ',iclm  c     if(distmin.le.clinc)then          !QUIQUI          
3560       $                 ,'( chi^2, ',RCHI2_STORE(ibest)                    
3561       $                 ,', norm.dist.= ', distmin                    CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<    
3562       $                 ,', cut ',clinc,' )'  *     ----------------------------
3563    c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
3564                      if(mod(VIEW(iclm),2).eq.0)then
3565                         XGOOD(nplanes-ip+1)=1.
3566                         resx(nplanes-ip+1)=rxmm
3567                         if(DEBUG.EQ.1)print*,'%%%% included X-cl ',iclm
3568         $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3569         $                    ,', dist.= ',distmin
3570         $                    ,', cut ',clinc,' )'
3571                      else
3572                         YGOOD(nplanes-ip+1)=1.
3573                         resy(nplanes-ip+1)=rymm
3574                         if(DEBUG.EQ.1)print*,'%%%% included Y-cl ',iclm
3575         $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3576         $                    ,', dist.= ', distmin
3577         $                    ,', cut ',clinc,' )'
3578                      endif
3579    c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
3580    *     ----------------------------
3581                      xm_A(nplanes-ip+1) = xmm_A
3582                      ym_A(nplanes-ip+1) = ymm_A
3583                      xm_B(nplanes-ip+1) = xmm_B
3584                      ym_B(nplanes-ip+1) = ymm_B
3585                      zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.
3586                      dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3587                      dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
3588    *     ----------------------------
3589                 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)  
 *              ----------------------------  
3590              endif              endif
3591           endif           endif
3592   133     continue   133     continue
# Line 2962  c              dedxtrk(nplanes-ip+1) = d Line 3605  c              dedxtrk(nplanes-ip+1) = d
3605  *                                                 *  *                                                 *
3606  *                                                 *  *                                                 *
3607  **************************************************  **************************************************
 cccccc 12/08/2006 modified by elena ---> (1)  
3608  *  *
3609        subroutine clean_XYclouds(ibest,iflag)        subroutine clean_XYclouds(ibest,iflag)
3610    
3611        include 'commontracker.f'        include 'commontracker.f'
3612        include 'level1.f'        include 'level1.f'
3613        include 'common_momanhough.f'        include 'common_momanhough.f'
3614  c      include 'momanhough_init.f'        include 'level2.f'      
       include 'level2.f'        !(1)  
 c      include 'calib.f'  
 c      include 'level1.f'  
   
3615    
3616          if(DEBUG.EQ.1)print*,'clean_XYclouds:'
3617    
3618        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3619    
# Line 2984  c      include 'level1.f' Line 3623  c      include 'level1.f'
3623              if(id.ne.0)then              if(id.ne.0)then
3624                 iclx=clx(ip,icp_cp(id))                 iclx=clx(ip,icp_cp(id))
3625                 icly=cly(ip,icp_cp(id))                 icly=cly(ip,icp_cp(id))
3626  c               cl_used(iclx)=1  !tag used clusters  c$$$               cl_used(iclx)=ntrk  !tag used clusters
3627  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)  
3628              elseif(icl.ne.0)then              elseif(icl.ne.0)then
3629  c               cl_used(icl)=1   !tag used clusters  c$$$               cl_used(icl)=ntrk   !tag used clusters
                cl_used(icl)=ntrk   !tag used clusters !1)  
3630              endif              endif
3631                            
 c               if(DEBUG)then  
 c                  print*,ip,' <<< ',id  
 c               endif  
3632  *     -----------------------------  *     -----------------------------
3633  *     remove the couple from clouds  *     remove the couple from clouds
3634  *     remove also vitual couples containing the  *     remove also vitual couples containing the
# Line 3012  c               endif Line 3645  c               endif
3645       $              cly(ip,icp).eq.icl       $              cly(ip,icp).eq.icl
3646       $              )then       $              )then
3647                    id=id_cp(ip,icp,1)                    id=id_cp(ip,icp,1)
3648                    if(DEBUG)then                    if(DEBUG.EQ.1)then
3649                       print*,ip,' <<< cp ',id                       print*,ip,' <<< cp ',id
3650       $                    ,' ( cl-x '       $                    ,' ( cl-x '
3651       $                    ,clx(ip,icp)       $                    ,clx(ip,icp)
# Line 3056  c               endif Line 3689  c               endif
3689        include 'level1.f'        include 'level1.f'
3690        include 'common_momanhough.f'        include 'common_momanhough.f'
3691        include 'level2.f'        include 'level2.f'
 c      include 'level1.f'  
3692    
3693    *     ---------------------------------
3694    *     variables initialized from level1
3695    *     ---------------------------------
3696        do i=1,nviews        do i=1,nviews
3697           good2(i)=good1(i)           good2(i)=good1(i)
3698             do j=1,nva1_view
3699                vkflag(i,j)=1
3700                if(cnnev(i,j).le.0)then
3701                   vkflag(i,j)=cnnev(i,j)
3702                endif
3703             enddo
3704        enddo        enddo
3705    *     ----------------
3706    *     level2 variables
3707    *     ----------------
3708        NTRK = 0        NTRK = 0
3709        do it=1,NTRKMAX        do it=1,NTRKMAX
3710           IMAGE(IT)=0           IMAGE(IT)=0
# Line 3073  c      include 'level1.f' Line 3715  c      include 'level1.f'
3715              ZM_nt(IP,IT) = 0              ZM_nt(IP,IT) = 0
3716              RESX_nt(IP,IT) = 0              RESX_nt(IP,IT) = 0
3717              RESY_nt(IP,IT) = 0              RESY_nt(IP,IT) = 0
3718                TAILX_nt(IP,IT) = 0
3719                TAILY_nt(IP,IT) = 0
3720                XBAD(IP,IT) = 0
3721                YBAD(IP,IT) = 0
3722              XGOOD_nt(IP,IT) = 0              XGOOD_nt(IP,IT) = 0
3723              YGOOD_nt(IP,IT) = 0              YGOOD_nt(IP,IT) = 0
3724                LS(IP,IT) = 0
3725              DEDX_X(IP,IT) = 0              DEDX_X(IP,IT) = 0
3726              DEDX_Y(IP,IT) = 0              DEDX_Y(IP,IT) = 0
3727              CLTRX(IP,IT) = 0              CLTRX(IP,IT) = 0
# Line 3207  c      include 'level1.f' Line 3854  c      include 'level1.f'
3854    
3855            
3856        include 'commontracker.f'        include 'commontracker.f'
 c      include 'level1.f'  
3857        include 'level1.f'        include 'level1.f'
3858        include 'common_momanhough.f'        include 'common_momanhough.f'
3859        include 'level2.f'        include 'level2.f'
3860        include 'common_mini_2.f'        include 'common_mini_2.f'
3861        real sinth,phi,pig              include 'calib.f'
3862    
3863          character*10 PFA
3864          common/FINALPFA/PFA
3865    
3866          real sinth,phi,pig
3867          integer ssensor,sladder
3868        pig=acos(-1.)        pig=acos(-1.)
3869    
3870    *     -------------------------------------
3871        chi2_nt(ntr)        = sngl(chi2)        chi2_nt(ntr)        = sngl(chi2)
3872        nstep_nt(ntr)       = nstep        nstep_nt(ntr)       = nstep
3873    *     -------------------------------------
3874        phi   = al(4)                  phi   = al(4)          
3875        sinth = al(3)                    sinth = al(3)            
3876        if(sinth.lt.0)then              if(sinth.lt.0)then      
# Line 3230  c      include 'level1.f' Line 3883  c      include 'level1.f'
3883       $     phi = phi + 2*pig         $     phi = phi + 2*pig  
3884        al(4) = phi                      al(4) = phi              
3885        al(3) = sinth                    al(3) = sinth            
   
3886        do i=1,5        do i=1,5
3887           al_nt(i,ntr)     = sngl(al(i))           al_nt(i,ntr)     = sngl(al(i))
3888           do j=1,5           do j=1,5
3889              coval(i,j,ntr) = sngl(cov(i,j))              coval(i,j,ntr) = sngl(cov(i,j))
3890           enddo           enddo
3891        enddo        enddo
3892          *     -------------------------------------      
3893        do ip=1,nplanes           ! loop on planes        do ip=1,nplanes           ! loop on planes
3894           xgood_nt(ip,ntr) = int(xgood(ip))           xgood_nt(ip,ntr) = int(xgood(ip))
3895           ygood_nt(ip,ntr) = int(ygood(ip))           ygood_nt(ip,ntr) = int(ygood(ip))
# Line 3246  c      include 'level1.f' Line 3898  c      include 'level1.f'
3898           zm_nt(ip,ntr)    = sngl(zm(ip))           zm_nt(ip,ntr)    = sngl(zm(ip))
3899           RESX_nt(IP,ntr)  = sngl(resx(ip))           RESX_nt(IP,ntr)  = sngl(resx(ip))
3900           RESY_nt(IP,ntr)  = sngl(resy(ip))           RESY_nt(IP,ntr)  = sngl(resy(ip))
3901             TAILX_nt(IP,ntr) = 0.
3902             TAILY_nt(IP,ntr) = 0.
3903           xv_nt(ip,ntr)    = sngl(xv(ip))           xv_nt(ip,ntr)    = sngl(xv(ip))
3904           yv_nt(ip,ntr)    = sngl(yv(ip))           yv_nt(ip,ntr)    = sngl(yv(ip))
3905           zv_nt(ip,ntr)    = sngl(zv(ip))           zv_nt(ip,ntr)    = sngl(zv(ip))
3906           axv_nt(ip,ntr)   = sngl(axv(ip))           axv_nt(ip,ntr)   = sngl(axv(ip))
3907           ayv_nt(ip,ntr)   = sngl(ayv(ip))           ayv_nt(ip,ntr)   = sngl(ayv(ip))  
3908           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)  
3909           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)             factor = sqrt(
3910         $        tan( acos(-1.) * sngl(axv(ip)) /180. )**2 +
3911         $        tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 +
3912         $        1. )
3913    
3914             dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)
3915             dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)  
3916        
3917           id  = CP_STORE(ip,IDCAND)           ax   = axv_nt(ip,ntr)
3918             ay   = ayv_nt(ip,ntr)
3919             bfx  = BX_STORE(ip,IDCAND)
3920             bfy  = BY_STORE(ip,IDCAND)
3921             if(ip.eq.6) ax = -1. * axv_nt(ip,ntr)
3922             if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND)
3923             tgtemp   = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001
3924             angx     = 180.*atan(tgtemp)/acos(-1.)
3925             tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
3926             angy    = 180.*atan(tgtemp)/acos(-1.)
3927            
3928    c         print*,'* ',ip,bfx,bfy,angx,angy
3929    
3930             id  = CP_STORE(ip,IDCAND) ! couple id
3931           icl = CLS_STORE(ip,IDCAND)           icl = CLS_STORE(ip,IDCAND)
3932             ssensor = -1
3933             sladder = -1
3934             ssensor = SENSOR_STORE(ip,IDCAND)
3935             sladder = LADDER_STORE(ip,IDCAND)
3936             if(ip.eq.6.and.ssensor.ne.0)ssensor = 3 - ssensor !notazione paolo x align
3937             LS(IP,ntr)      = ssensor+10*sladder
3938    
3939           if(id.ne.0)then           if(id.ne.0)then
3940    c           >>> is a couple
3941              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))
3942              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))
3943  c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)  
3944                cl_used(cltrx(ip,ntr)) = 1     !tag used clusters          
3945                cl_used(cltry(ip,ntr)) = 1     !tag used clusters          
3946    
3947    c$$$            nnnnx = npfastrips(clx(nplanes-ip+1,icp_cp(id)),PFA,angx)
3948    c$$$            nnnny = npfastrips(cly(nplanes-ip+1,icp_cp(id)),PFA,angy)            
3949    c$$$            xbad(ip,ntr)= nbadstrips(nnnnx,clx(nplanes-ip+1,icp_cp(id)))
3950    c$$$            ybad(ip,ntr)= nbadstrips(nnnny,cly(nplanes-ip+1,icp_cp(id)))
3951                xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))
3952                ybad(ip,ntr)= nbadstrips(4,cly(nplanes-ip+1,icp_cp(id)))
3953    
3954    
3955                if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0)
3956         $           dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3957                if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)
3958         $           dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3959    
3960           elseif(icl.ne.0)then           elseif(icl.ne.0)then
3961              if(mod(VIEW(icl),2).eq.0)cltrx(ip,ntr)=icl  
3962              if(mod(VIEW(icl),2).eq.1)cltry(ip,ntr)=icl              cl_used(icl) = 1    !tag used clusters          
3963  c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)  
3964                if(mod(VIEW(icl),2).eq.0)then
3965                   cltrx(ip,ntr)=icl
3966    c$$$               nnnnn = npfastrips(icl,PFA,angx)
3967    c$$$               xbad(ip,ntr) = nbadstrips(nnnnn,icl)
3968                   xbad(ip,ntr) = nbadstrips(4,icl)
3969    
3970                   if(nsatstrips(icl).gt.0)dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3971                elseif(mod(VIEW(icl),2).eq.1)then
3972                   cltry(ip,ntr)=icl
3973    c$$$               nnnnn = npfastrips(icl,PFA,angy)
3974    c$$$               ybad(ip,ntr) = nbadstrips(nnnnn,icl)
3975                   ybad(ip,ntr) = nbadstrips(4,icl)
3976                   if(nsatstrips(icl).gt.0)dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3977                endif
3978    
3979           endif                     endif          
3980    
3981        enddo        enddo
3982    
3983          if(DEBUG.eq.1)then
3984             print*,'> STORING TRACK ',ntr
3985             print*,'clusters: '
3986             do ip=1,6
3987                print*,'> ',ip,' -- ',cltrx(ip,ntr),cltry(ip,ntr)
3988             enddo
3989          endif
3990    
3991    c$$$      print*,(xgood(i),i=1,6)
3992    c$$$      print*,(ygood(i),i=1,6)
3993    c$$$      print*,(ls(i,ntr),i=1,6)
3994    c$$$      print*,(dedx_x(i,ntr),i=1,6)
3995    c$$$      print*,(dedx_y(i,ntr),i=1,6)
3996    c$$$      print*,'-----------------------'
3997    
3998        end        end
3999    
# Line 3280  c            print*,ip,' ',cltrx(ip,ntr) Line 4006  c            print*,ip,' ',cltrx(ip,ntr)
4006  *     -------------------------------------------------------  *     -------------------------------------------------------
4007    
4008        include 'commontracker.f'        include 'commontracker.f'
 c      include 'level1.f'  
4009        include 'calib.f'        include 'calib.f'
4010        include 'level1.f'        include 'level1.f'
4011        include 'common_momanhough.f'        include 'common_momanhough.f'
# Line 3288  c      include 'level1.f' Line 4013  c      include 'level1.f'
4013        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
4014    
4015  *     count #cluster per plane not associated to any track  *     count #cluster per plane not associated to any track
 c      good2=1!.true.  
4016        nclsx = 0        nclsx = 0
4017        nclsy = 0        nclsy = 0
4018    
4019        do iv = 1,nviews        do iv = 1,nviews
4020           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)
4021             good2(iv) = good2(iv) + mask_view(iv)*2**8
4022        enddo        enddo
4023    
4024          if(DEBUG.eq.1)then
4025             print*,'> STORING SINGLETS '
4026          endif
4027    
4028        do icl=1,nclstr1        do icl=1,nclstr1
4029    
4030             ip=nplanes-npl(VIEW(icl))+1            
4031            
4032           if(cl_used(icl).eq.0)then !cluster not included in any track           if(cl_used(icl).eq.0)then !cluster not included in any track
             ip=nplanes-npl(VIEW(icl))+1              
4033              if(mod(VIEW(icl),2).eq.0)then !=== X views              if(mod(VIEW(icl),2).eq.0)then !=== X views
4034                 nclsx = nclsx + 1                 nclsx = nclsx + 1
4035                 planex(nclsx) = ip                 planex(nclsx) = ip
4036                 sgnlxs(nclsx) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
4037                   if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx)
4038                 clsx(nclsx)   = icl                 clsx(nclsx)   = icl
4039                 do is=1,2                 do is=1,2
4040  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
4041                    call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)  c                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
4042                      call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.,0.,0.)
4043                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2
4044                 enddo                 enddo
4045  c$$$               print*,'nclsx         ',nclsx  c$$$               print*,'nclsx         ',nclsx
# Line 3317  c$$$               print*,'xs(2,nclsx)   Line 4050  c$$$               print*,'xs(2,nclsx)  
4050              else                          !=== Y views              else                          !=== Y views
4051                 nclsy = nclsy + 1                 nclsy = nclsy + 1
4052                 planey(nclsy) = ip                 planey(nclsy) = ip
4053                 sgnlys(nclsy) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
4054                   if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy)
4055                 clsy(nclsy)   = icl                 clsy(nclsy)   = icl
4056                 do is=1,2                 do is=1,2
4057  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
4058                    call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)  c                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
4059                      call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.,0.,0.)
4060                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2
4061                 enddo                 enddo
4062  c$$$               print*,'nclsy         ',nclsy  c$$$               print*,'nclsy         ',nclsy
# Line 3331  c$$$               print*,'ys(1,nclsy)   Line 4066  c$$$               print*,'ys(1,nclsy)  
4066  c$$$               print*,'ys(2,nclsy)   ',ys(2,nclsy)  c$$$               print*,'ys(2,nclsy)   ',ys(2,nclsy)
4067              endif              endif
4068           endif           endif
 c      print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip)  
4069    
4070  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO
4071           whichtrack(icl) = cl_used(icl)           whichtrack(icl) = cl_used(icl)
4072    *     --------------------------------------------------
4073    *     per non perdere la testa...
4074    *     whichtrack(icl) e` una variabile del common level1
4075    *     che serve solo per sapere quali cluster sono stati
4076    *     associati ad una traccia, e permettere di salvare
4077    *     solo questi nell'albero di uscita
4078    *     --------------------------------------------------
4079            
4080    
4081    c$$$         print*,' cl ',icl,' --> ',cl_used(icl)
4082    c$$$
4083    c$$$         if( cl_used(icl).ne.0 )then
4084    c$$$            if(
4085    c$$$     $           mod(VIEW(icl),2).eq.0.and.
4086    c$$$     $           cltrx(ip,whichtrack(icl)).ne.icl )
4087    c$$$     $           print*,'**WARNING** cltrx(',ip,',',whichtrack(icl)
4088    c$$$     $           ,')=',cltrx(ip,whichtrack(icl)),'.ne.',icl
4089    c$$$            if(
4090    c$$$     $           mod(VIEW(icl),2).eq.1.and.
4091    c$$$     $           cltry(ip,whichtrack(icl)).ne.icl )
4092    c$$$     $           print*,'**WARNING** cltry(',ip,',',whichtrack(icl)
4093    c$$$     $           ,')=',cltry(ip,whichtrack(icl)),'.ne.',icl
4094    c$$$         endif
4095            
4096    
4097        enddo        enddo
4098        end        end
Legend:
Removed from v.1.16  
changed lines
  Added in v.1.28
  ViewVC Help
Powered by ViewVC 1.1.23