/[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.13 by pam-fi, Wed Nov 8 16:42:28 2006 UTC revision 1.19 by pam-fi, Tue Feb 20 08:32:52 2007 UTC
# Line 18  Line 18 
18        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
19        include 'common_mini_2.f'        include 'common_mini_2.f'
20        include 'calib.f'        include 'calib.f'
 c      include 'level1.f'  
21        include 'level2.f'        include 'level2.f'
22    
23  c      include 'momanhough_init.f'  c      include 'momanhough_init.f'
# Line 234  c$$$            endif Line 233  c$$$            endif
233  c$$$         enddo  c$$$         enddo
234  c$$$         if(ibest.eq.0)goto 880 !>> no good candidates  c$$$         if(ibest.eq.0)goto 880 !>> no good candidates
235    
236           rchi2best=1000000000.  *     -------------------------------------------------------
237    *     order track-candidates according to:
238    *     1st) decreasing n.points
239    *     2nd) increasing chi**2
240    *     -------------------------------------------------------
241             rchi2best=1000000000.
242           ndofbest=0             !(1)           ndofbest=0             !(1)
243           do i=1,ntracks           do i=1,ntracks
244             if(RCHI2_STORE(i).lt.rchi2best.and.             ndof=0               !(1)
245       $          RCHI2_STORE(i).gt.0)then             do ii=1,nplanes      !(1)
246               ndof=0             !(1)               ndof=ndof          !(1)
247               do ii=1,nplanes    !(1)       $            +int(xgood_store(ii,i)) !(1)
248                 ndof=ndof        !(1)       $            +int(ygood_store(ii,i)) !(1)
249       $              +int(xgood_store(ii,i)) !(1)             enddo                !(1)
250       $              +int(ygood_store(ii,i)) !(1)             if(ndof.gt.ndofbest)then !(1)
251               enddo              !(1)               ibest=i
252               if(ndof.ge.ndofbest)then !(1)               rchi2best=RCHI2_STORE(i)
253                 ndofbest=ndof      !(1)
254               elseif(ndof.eq.ndofbest)then !(1)
255                 if(RCHI2_STORE(i).lt.rchi2best.and.
256         $            RCHI2_STORE(i).gt.0)then
257                 ibest=i                 ibest=i
258                 rchi2best=RCHI2_STORE(i)                 rchi2best=RCHI2_STORE(i)
259                 ndofbest=ndof    !(1)                 ndofbest=ndof    !(1)
260               endif              !(1)               endif              !(1)
261             endif             endif
262           enddo           enddo
263    
264    c$$$         rchi2best=1000000000.
265    c$$$         ndofbest=0             !(1)
266    c$$$         do i=1,ntracks
267    c$$$           if(RCHI2_STORE(i).lt.rchi2best.and.
268    c$$$     $          RCHI2_STORE(i).gt.0)then
269    c$$$             ndof=0             !(1)
270    c$$$             do ii=1,nplanes    !(1)
271    c$$$               ndof=ndof        !(1)
272    c$$$     $              +int(xgood_store(ii,i)) !(1)
273    c$$$     $              +int(ygood_store(ii,i)) !(1)
274    c$$$             enddo              !(1)
275    c$$$             if(ndof.ge.ndofbest)then !(1)
276    c$$$               ibest=i
277    c$$$               rchi2best=RCHI2_STORE(i)
278    c$$$               ndofbest=ndof    !(1)
279    c$$$             endif              !(1)
280    c$$$           endif
281    c$$$         enddo
282    
283           if(ibest.eq.0)goto 880 !>> no good candidates           if(ibest.eq.0)goto 880 !>> no good candidates
284  *-------------------------------------------------------------------------------      *-------------------------------------------------------------------------------    
285  *     The best track candidate (ibest) is selected and a new fitting is performed.  *     The best track candidate (ibest) is selected and a new fitting is performed.
# Line 286  c$$$         if(ibest.eq.0)goto 880 !>> Line 314  c$$$         if(ibest.eq.0)goto 880 !>>
314           do i=1,5           do i=1,5
315              AL_GUESS(i)=AL(i)              AL_GUESS(i)=AL(i)
316           enddo           enddo
317    c         print*,'## guess: ',al
318    
319           do i=1,5           do i=1,5
320              AL(i)=dble(AL_STORE(i,icand))                          AL(i)=dble(AL_STORE(i,icand))            
# Line 298  c$$$         if(ibest.eq.0)goto 880 !>> Line 327  c$$$         if(ibest.eq.0)goto 880 !>>
327           iprint=0           iprint=0
328  c         if(DEBUG)iprint=1  c         if(DEBUG)iprint=1
329           if(VERBOSE)iprint=1           if(VERBOSE)iprint=1
330             if(DEBUG)iprint=2
331           call mini2(jstep,ifail,iprint)           call mini2(jstep,ifail,iprint)
332           if(ifail.ne.0) then           if(ifail.ne.0) then
333              if(.true.)then              if(VERBOSE)then
334                 print *,                 print *,
335       $              '*** MINIMIZATION FAILURE *** (after refinement) '       $              '*** MINIMIZATION FAILURE *** (after refinement) '
336       $              ,iev       $              ,iev
# Line 423  c     $        rchi2best.lt.15..and. Line 453  c     $        rchi2best.lt.15..and.
453        end        end
454    
455    
   
   
 c$$$************************************************************  
 c$$$  
 c$$$      subroutine readmipparam  
 c$$$              
 c$$$      include 'commontracker.f'  
 c$$$      include 'calib.f'  
 c$$$  
 c$$$      character*60 fname_param  
 c$$$ 201  format('trk-LADDER',i1,'-mip.dat')  
 c$$$      do ilad=1,nladders_view          
 c$$$         write(fname_param,201)ilad  
 c$$$         print *,'Opening file: ',fname_param  
 c$$$         open(10,  
 c$$$     $        FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param))  
 c$$$     $        ,STATUS='UNKNOWN'  
 c$$$     $        ,IOSTAT=iostat  
 c$$$     $        )  
 c$$$         if(iostat.ne.0)then  
 c$$$            print*,'READMIPPARAM: *** Error in opening file ***'  
 c$$$            return  
 c$$$         endif  
 c$$$         do iv=1,nviews  
 c$$$            read(10,*  
 c$$$     $           ,IOSTAT=iostat  
 c$$$     $           )pip,  
 c$$$     $            mip(int(pip),ilad)  
 c$$$c            print*,ilad,iv,pip,mip(int(pip),ilad)  
 c$$$         enddo  
 c$$$         close(10)  
 c$$$      enddo  
 c$$$  
 c$$$      return  
 c$$$      end  
 c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  
 c$$$      subroutine readchargeparam  
 c$$$        
 c$$$        
 c$$$      include 'commontracker.f'  
 c$$$      include 'calib.f'  
 c$$$  
 c$$$      character*60 fname_param  
 c$$$ 201  format('charge-l',i1,'.dat')  
 c$$$      do ilad=1,nladders_view          
 c$$$         write(fname_param,201)ilad  
 c$$$         print *,'Opening file: ',fname_param  
 c$$$         open(10,  
 c$$$     $        FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param))  
 c$$$     $        ,STATUS='UNKNOWN'  
 c$$$     $        ,IOSTAT=iostat  
 c$$$     $        )  
 c$$$         if(iostat.ne.0)then  
 c$$$            print*,'READCHARGEPARAM: *** Error in opening file ***'  
 c$$$            return  
 c$$$         endif  
 c$$$         do ip=1,nplanes  
 c$$$            read(10,*  
 c$$$     $           ,IOSTAT=iostat  
 c$$$     $           )pip,  
 c$$$     $            kch(ip,ilad),cch(ip,ilad),sch(ip,ilad)          
 c$$$c            print*,ilad,ip,pip,kch(ip,ilad),  
 c$$$c     $           cch(ip,ilad),sch(ip,ilad)  
 c$$$         enddo  
 c$$$         close(10)  
 c$$$      enddo  
 c$$$  
 c$$$      return  
 c$$$      end  
 c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  
 c$$$      subroutine readetaparam  
 c$$$*     -----------------------------------------  
 c$$$*     read eta2,3,4 calibration parameters  
 c$$$*     and fill variables:  
 c$$$*  
 c$$$*     eta2(netabin,nladders_view,nviews)  
 c$$$*     eta3(2*netabin,nladders_view,nviews)  
 c$$$*     eta4(2*netabin,nladders_view,nviews)  
 c$$$*  
 c$$$      include 'commontracker.f'  
 c$$$      include 'calib.f'  
 c$$$  
 c$$$      character*40 fname_binning  
 c$$$      character*40 fname_param  
 c$$$c      character*120 cmd1  
 c$$$c      character*120 cmd2  
 c$$$  
 c$$$  
 c$$$******retrieve ANGULAR BINNING info  
 c$$$      fname_binning='binning.dat'  
 c$$$      print *,'Opening file: ',fname_binning  
 c$$$      open(10,  
 c$$$     $     FILE='./bin-aux/'//fname_binning(1:LNBLNK(fname_binning))  
 c$$$     $     ,STATUS='UNKNOWN'  
 c$$$     $     ,IOSTAT=iostat  
 c$$$     $     )  
 c$$$      if(iostat.ne.0)then  
 c$$$         print*,'READETAPARAM: *** Error in opening file ***'  
 c$$$         return  
 c$$$      endif  
 c$$$      print*,'---- ANGULAR BINNING ----'  
 c$$$      print*,'Bin   -   angL   -   angR'  
 c$$$ 101  format(i2,'       ',f6.2,'     ',f6.2)  
 c$$$      do ibin=1,nangmax  
 c$$$         read(10,*  
 c$$$     $        ,IOSTAT=iostat  
 c$$$     $        )xnn,angL(ibin),angR(ibin)  
 c$$$         if(iostat.ne.0)goto 1000  
 c$$$         write(*,101)int(xnn),angL(ibin),angR(ibin)  
 c$$$      enddo          
 c$$$ 1000 nangbin=int(xnn)  
 c$$$      close(10)  
 c$$$      print*,'-------------------------'  
 c$$$        
 c$$$  
 c$$$  
 c$$$      do ieta=2,4               !loop on eta 2,3,4          
 c$$$******retrieve correction parameters  
 c$$$ 200     format(' Opening eta',i1,' files...')  
 c$$$         write(*,200)ieta  
 c$$$  
 c$$$ 201     format('eta',i1,'-bin',i1,'-l',i1,'.dat')  
 c$$$ 202     format('eta',i1,'-bin',i2,'-l',i1,'.dat')  
 c$$$         do iang=1,nangbin  
 c$$$            do ilad=1,nladders_view  
 c$$$               if(iang.lt.10)write(fname_param,201)ieta,iang,ilad  
 c$$$               if(iang.ge.10)write(fname_param,202)ieta,iang,ilad  
 c$$$c               print *,'Opening file: ',fname_param  
 c$$$               open(10,  
 c$$$     $             FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param))  
 c$$$     $              ,STATUS='UNKNOWN'  
 c$$$     $              ,IOSTAT=iostat  
 c$$$     $              )  
 c$$$               if(iostat.ne.0)then  
 c$$$                  print*,'READETAPARAM: *** Error in opening file ***'  
 c$$$                  return  
 c$$$               endif  
 c$$$               do ival=1,netavalmax  
 c$$$                  if(ieta.eq.2)read(10,*  
 c$$$     $                 ,IOSTAT=iostat  
 c$$$     $                 )  
 c$$$     $                 eta2(ival,iang),  
 c$$$     $                 (feta2(ival,iv,ilad,iang),iv=1,nviews)  
 c$$$                  if(ieta.eq.3)read(10,*  
 c$$$     $                 ,IOSTAT=iostat  
 c$$$     $                 )  
 c$$$     $                 eta3(ival,iang),  
 c$$$     $                 (feta3(ival,iv,ilad,iang),iv=1,nviews)  
 c$$$                  if(ieta.eq.4)read(10,*  
 c$$$     $                 ,IOSTAT=iostat  
 c$$$     $                 )  
 c$$$     $                 eta4(ival,iang),  
 c$$$     $                 (feta4(ival,iv,ilad,iang),iv=1,nviews)  
 c$$$                  if(iostat.ne.0)then  
 c$$$                     netaval=ival-1  
 c$$$c$$$                     if(eta2(1,iang).ne.-eta2(netaval,iang))  
 c$$$c$$$     $                    print*,'**** ERROR on parameters !!! ****'  
 c$$$                     goto 2000  
 c$$$                  endif  
 c$$$               enddo  
 c$$$ 2000          close(10)  
 c$$$*               print*,'... done'  
 c$$$            enddo  
 c$$$         enddo  
 c$$$  
 c$$$      enddo                     !end loop on eta 2,3,4  
 c$$$  
 c$$$  
 c$$$      return  
 c$$$      end  
 c$$$  
   
456                
457  ************************************************************  ************************************************************
458  ************************************************************  ************************************************************
# Line 619  c$$$ Line 477  c$$$
477  *     PFAy   - Position Finding Algorithm in y (COG2,ETA2,...)  *     PFAy   - Position Finding Algorithm in y (COG2,ETA2,...)
478  *     angx   - Projected angle in x  *     angx   - Projected angle in x
479  *     angy   - Projected angle in y  *     angy   - Projected angle in y
480    *     bfx    - x-component of magnetci field
481    *     bfy    - y-component of magnetic field
482  *  *
483  *     --------- COUPLES -------------------------------------------------------  *     --------- COUPLES -------------------------------------------------------
484  *     The couple defines a point in the space.  *     The couple defines a point in the space.
# Line 657  c$$$ Line 517  c$$$
517  *  *
518  *  *
519    
520        subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,angx,angy)        subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
521    
 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*****************************************************  
522                
523        include 'commontracker.f'        include 'commontracker.f'
524        include 'level1.f'        include 'level1.f'
# Line 684  c      common/dbg/DEBUG Line 536  c      common/dbg/DEBUG
536        integer sensor        integer sensor
537        integer viewx,viewy        integer viewx,viewy
538        character*4 PFAx,PFAy     !PFA to be used        character*4 PFAx,PFAy     !PFA to be used
539        real angx,angy            !X-Y angle        real ax,ay                !X-Y geometric angle
540          real angx,angy            !X-Y effective angle
541          real bfx,bfy              !X-Y b-field components
542    
543        real stripx,stripy        real stripx,stripy
544    
# Line 710  c      double precision xi_B,yi_B,zi_B Line 564  c      double precision xi_B,yi_B,zi_B
564        xPAM_B = 0.        xPAM_B = 0.
565        yPAM_B = 0.        yPAM_B = 0.
566        zPAM_B = 0.        zPAM_B = 0.
567    c      print*,'## xyz_PAM: ',icx,icy,sensor,PFAx,PFAy,angx,angy
568  *     -----------------  *     -----------------
569  *     CLUSTER X  *     CLUSTER X
570  *     -----------------  *     -----------------
571    
572        if(icx.ne.0)then        if(icx.ne.0)then
573    
574           viewx = VIEW(icx)           viewx = VIEW(icx)
575           nldx = nld(MAXS(icx),VIEW(icx))           nldx = nld(MAXS(icx),VIEW(icx))
576           nplx = npl(VIEW(icx))           nplx = npl(VIEW(icx))
577           resxPAM = RESXAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!           resxPAM = RESXAV
           
578           stripx = float(MAXS(icx))           stripx = float(MAXS(icx))
579           if(PFAx.eq.'COG1')then  !(1)  *        --------------------------
580              stripx = stripx      !(1)  *        magnetic-field corrections
581              resxPAM = resxPAM    !(1)  *        --------------------------
582    c$$$         print*,nplx,ax,bfy/10.
583             angtemp  = ax
584             bfytemp  = bfy
585             if(nplx.eq.6) angtemp = -1. * ax
586             if(nplx.eq.6) bfytemp = -1. * bfy
587             tgtemp = tan(angtemp*acos(-1.)/180.) + pmuH_h*bfytemp*0.00001
588             angx = 180.*atan(tgtemp)/acos(-1.)
589             stripx = stripx - 0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX
590    c$$$         print*,angx,0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX
591    c$$$         print*,'========================'
592    *        --------------------------
593    
594             if(PFAx.eq.'COG1')then
595                stripx = stripx    
596                resxPAM = resxPAM  
597           elseif(PFAx.eq.'COG2')then           elseif(PFAx.eq.'COG2')then
598              stripx = stripx + cog(2,icx)                          stripx = stripx + cog(2,icx)            
599              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
600             elseif(PFAx.eq.'COG3')then
601                stripx = stripx + cog(3,icx)            
602                resxPAM = resxPAM*fbad_cog(3,icx)
603             elseif(PFAx.eq.'COG4')then
604    c            print*,'COG4'
605                stripx = stripx + cog(4,icx)            
606                resxPAM = resxPAM*fbad_cog(4,icx)
607           elseif(PFAx.eq.'ETA2')then           elseif(PFAx.eq.'ETA2')then
608  c            cog2 = cog(2,icx)              stripx = stripx + pfaeta2(icx,angx)          
609  c            etacorr = pfaeta2(cog2,viewx,nldx,angx)                          resxPAM = risx_eta2(abs(angx))                    
 c            stripx = stripx + etacorr  
             stripx = stripx + pfaeta2(icx,angx)            !(3)  
             resxPAM = risx_eta2(angx)                       !   (4)  
610              if(DEBUG.and.fbad_cog(2,icx).ne.1)              if(DEBUG.and.fbad_cog(2,icx).ne.1)
611       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
612              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
613           elseif(PFAx.eq.'ETA3')then                         !(3)           elseif(PFAx.eq.'ETA3')then                        
614              stripx = stripx + pfaeta3(icx,angx)            !(3)              stripx = stripx + pfaeta3(icx,angx)          
615              resxPAM = risx_eta3(angx)                       !   (4)              resxPAM = risx_eta3(abs(angx))                      
616              if(DEBUG.and.fbad_cog(3,icx).ne.1)              !(3)              if(DEBUG.and.fbad_cog(3,icx).ne.1)            
617       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)
618              resxPAM = resxPAM*fbad_cog(3,icx)               !(3)              resxPAM = resxPAM*fbad_cog(3,icx)              
619           elseif(PFAx.eq.'ETA4')then                         !(3)           elseif(PFAx.eq.'ETA4')then                        
620              stripx = stripx + pfaeta4(icx,angx)            !(3)              stripx = stripx + pfaeta4(icx,angx)            
621              resxPAM = risx_eta4(angx)                       !   (4)              resxPAM = risx_eta4(abs(angx))                      
622              if(DEBUG.and.fbad_cog(4,icx).ne.1)              !(3)              if(DEBUG.and.fbad_cog(4,icx).ne.1)              
623       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)
624              resxPAM = resxPAM*fbad_cog(4,icx)               !(3)              resxPAM = resxPAM*fbad_cog(4,icx)              
625           elseif(PFAx.eq.'ETA')then                          !(3)           elseif(PFAx.eq.'ETA')then  
626              stripx = stripx + pfaeta(icx,angx)             !(3)  c            print*,'ETA'
627              resxPAM = ris_eta(icx,angx)                     !   (4)              stripx = stripx + pfaeta(icx,angx)            
628              if(DEBUG.and.fbad_cog(2,icx).ne.1)              !(3)              resxPAM = ris_eta(icx,angx)                    
629       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)              if(DEBUG.and.fbad_cog(2,icx).ne.1)              
630  c            resxPAM = resxPAM*fbad_cog(2,icx)               !(3)TEMPORANEO       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
631              resxPAM = resxPAM*fbad_eta(icx,angx)             !(3)(4)              resxPAM = resxPAM*fbad_eta(icx,angx)            
632           elseif(PFAx.eq.'COG')then           !(2)           elseif(PFAx.eq.'COG')then          
633              stripx = stripx + cog(0,icx)     !(2)                      stripx = stripx + cog(0,icx)            
634              resxPAM = risx_cog(angx)                        !   (4)              resxPAM = risx_cog(abs(angx))                    
635              resxPAM = resxPAM*fbad_cog(0,icx)!(2)              resxPAM = resxPAM*fbad_cog(0,icx)
636           else           else
637              print*,'*** Non valid p.f.a. (x) --> ',PFAx              print*,'*** Non valid p.f.a. (x) --> ',PFAx
638           endif           endif
639    
640    c         print*,'%%%%%%%%%%%%'
641    
642        endif        endif
 c      if(icy.eq.0.and.icx.ne.0)  
 c     $     print*,PFAx,icx,angx,stripx,resxPAM,'***'  
643                
644  *     -----------------  *     -----------------
645  *     CLUSTER Y  *     CLUSTER Y
646  *     -----------------  *     -----------------
647    
648        if(icy.ne.0)then        if(icy.ne.0)then
649    
650           viewy = VIEW(icy)           viewy = VIEW(icy)
651           nldy = nld(MAXS(icy),VIEW(icy))           nldy = nld(MAXS(icy),VIEW(icy))
652           nply = npl(VIEW(icy))           nply = npl(VIEW(icy))
653           resyPAM = RESYAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!           resyPAM = RESYAV
654             stripy = float(MAXS(icy))
655    
656           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
657              print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '              print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '
658       $           ,icx,icy       $           ,icx,icy
659              goto 100              goto 100
660           endif           endif
661    *        --------------------------
662    *        magnetic-field corrections
663    *        --------------------------
664             tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
665             angy    = 180.*atan(tgtemp)/acos(-1.)
666             stripy = stripy + 0.5*pmuH_e*bfx*0.00001*SiDimZ/pitchY
667    *        --------------------------
668                    
          stripy = float(MAXS(icy))  
669           if(PFAy.eq.'COG1')then !(1)           if(PFAy.eq.'COG1')then !(1)
670              stripy = stripy     !(1)              stripy = stripy     !(1)
671              resyPAM = resyPAM   !(1)              resyPAM = resyPAM   !(1)
672           elseif(PFAy.eq.'COG2')then           elseif(PFAy.eq.'COG2')then
673              stripy = stripy + cog(2,icy)              stripy = stripy + cog(2,icy)
674              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
675             elseif(PFAy.eq.'COG3')then
676                stripy = stripy + cog(3,icy)
677                resyPAM = resyPAM*fbad_cog(3,icy)
678             elseif(PFAy.eq.'COG4')then
679                stripy = stripy + cog(4,icy)
680                resyPAM = resyPAM*fbad_cog(4,icy)
681           elseif(PFAy.eq.'ETA2')then           elseif(PFAy.eq.'ETA2')then
682  c            cog2 = cog(2,icy)  c            cog2 = cog(2,icy)
683  c            etacorr = pfaeta2(cog2,viewy,nldy,angy)  c            etacorr = pfaeta2(cog2,viewy,nldy,angy)
684  c            stripy = stripy + etacorr  c            stripy = stripy + etacorr
685              stripy = stripy + pfaeta2(icy,angy)            !(3)              stripy = stripy + pfaeta2(icy,angy)            !(3)
686              resyPAM = risy_eta2(angy)                       !   (4)              resyPAM = risy_eta2(abs(angy))                       !   (4)
687              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
688              if(DEBUG.and.fbad_cog(2,icy).ne.1)              if(DEBUG.and.fbad_cog(2,icy).ne.1)
689       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
# Line 820  c            resyPAM = resyPAM*fbad_cog( Line 706  c            resyPAM = resyPAM*fbad_cog(
706       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)!(3)       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)!(3)
707           elseif(PFAy.eq.'COG')then           elseif(PFAy.eq.'COG')then
708              stripy = stripy + cog(0,icy)                          stripy = stripy + cog(0,icy)            
709              resyPAM = risy_cog(angy)                        !   (4)              resyPAM = risy_cog(abs(angy))                        !   (4)
710  c            resyPAM = ris_eta(icy,angy)                    !   (4)  c            resyPAM = ris_eta(icy,angy)                    !   (4)
711              resyPAM = resyPAM*fbad_cog(0,icy)              resyPAM = resyPAM*fbad_cog(0,icy)
712           else           else
# Line 829  c            resyPAM = ris_eta(icy,angy) Line 715  c            resyPAM = ris_eta(icy,angy)
715    
716        endif        endif
717    
718          c      print*,'## stripx,stripy ',stripx,stripy
719    
720  c===========================================================  c===========================================================
721  C     COUPLE  C     COUPLE
722  C===========================================================  C===========================================================
# Line 1031  c         print*,'A-(',xPAM_A,yPAM_A,') Line 918  c         print*,'A-(',xPAM_A,yPAM_A,')
918                            
919        endif        endif
920                    
921    
922    c      print*,'## xPAM,yPAM,zPAM       ',xPAM,yPAM,zPAM
923    c      print*,'## xPAM_A,yPAM_A,zPAM_A ',xPAM_A,yPAM_A,zPAM_A
924    c      print*,'## xPAM_B,yPAM_B,zPAM_B ',xPAM_B,yPAM_B,zPAM_B
925    
926   100  continue   100  continue
927        end        end
928    
# Line 1444  c      include 'common_analysis.f' Line 1336  c      include 'common_analysis.f'
1336  *************************************************************************  *************************************************************************
1337  *************************************************************************  *************************************************************************
1338  *************************************************************************  *************************************************************************
 c$$$      subroutine book_debug  
 c$$$  
 c$$$      include 'commontracker.f'  
 c$$$      include 'common_momanhough.f'  
 c$$$      include 'common_level2debug.f'  
 c$$$  
 c$$$      character*35 block1,block2,block3!,block4  
 c$$$     $     ,block5!,block6        
 c$$$        
 c$$$* * * * * * * * * * * * * * * * * * * * * * * *  
 c$$$*     HOUGH TRANSFORM PARAMETERS  
 c$$$        
 c$$$      call HBOOK2(1003  
 c$$$     $     ,'y vs tg thyz'  
 c$$$     $     ,300,-1.,1.         !x          
 c$$$     $     ,3000,-70.,70.,0.)   !y  
 c$$$  
 c$$$      call HBOOK1(1004  
 c$$$     $     ,'Dy'  
 c$$$     $     ,100,0.,2.,0.)    
 c$$$  
 c$$$      call HBOOK1(1005  
 c$$$     $     ,'D thyz'  
 c$$$     $     ,100,0.,.05,0.)    
 c$$$  
 c$$$  
 c$$$  
 c$$$*     DEBUG ntuple:  
 c$$$      call HBNT(ntp_level2+1,'LEVEL2',' ')  
 c$$$        
 c$$$      call HBNAME(ntp_level2+1,'EVENT',good2_nt,  
 c$$$     $     'GOOD2:L,NEV2:I')  
 c$$$  
 c$$$ 411  format('NDBLT:I::[0,',I5,']')  
 c$$$      write(block1,411) ndblt_max_nt  
 c$$$      call HBNAME(ntp_level2+1,'HOUGH YZ',ndblt_nt,  
 c$$$     $     block1//'  
 c$$$     $     ,ALFAYZ1(NDBLT):R  
 c$$$     $     ,ALFAYZ2(NDBLT):R  
 c$$$     $     ,DB_CLOUD(NDBLT):I  
 c$$$     $     ')    
 c$$$  
 c$$$ 412  format('NTRPT:I::[0,',I5,']')  
 c$$$      write(block2,412) ntrpt_max_nt  
 c$$$      call HBNAME(ntp_level2+1,'HOUGH XZ',NTRPT_nt,  
 c$$$     $     block2//'  
 c$$$     $     ,ALFAXZ1(NTRPT):R  
 c$$$     $     ,ALFAXZ2(NTRPT):R    
 c$$$     $     ,ALFAXZ3(NTRPT):R  
 c$$$     $     ,TR_CLOUD(NTRPT):I  
 c$$$     $     ')  
 c$$$        
 c$$$    
 c$$$ 413  format('NCLOUDS_YZ:I::[0,',I4,']')  
 c$$$c$$$ 414  format('DB_CLOUD(',I4,'):I')  
 c$$$      write(block3,413) ncloyz_max  
 c$$$c$$$      write(block4,414) ndblt_max_nt  
 c$$$      call HBNAME(ntp_level2+1,'CLOUD YZ',NCLOUDS_YZ,  
 c$$$     $     block3//'  
 c$$$     $     ,ALFAYZ1_AV(NCLOUDS_YZ):R  
 c$$$     $     ,ALFAYZ2_AV(NCLOUDS_YZ):R  
 c$$$     $     ,PTCLOUD_YZ(NCLOUDS_YZ):I'  
 c$$$c$$$     $     ,'//block4  
 c$$$     $     )  
 c$$$  
 c$$$ 415  format('NCLOUDS_XZ:I::[0,',I4,']')  
 c$$$c$$$ 416  format('TR_CLOUD(',I5,'):I')  
 c$$$      write(block5,415) ncloxz_max  
 c$$$c$$$      write(block6,416) ntrpt_max_nt  
 c$$$      call HBNAME(ntp_level2+1,'CLOUD XZ',NCLOUDS_XZ,  
 c$$$     $     block5//'  
 c$$$     $     ,ALFAXZ1_AV(NCLOUDS_XZ):R  
 c$$$     $     ,ALFAXZ2_AV(NCLOUDS_XZ):R  
 c$$$     $     ,ALFAXZ3_AV(NCLOUDS_XZ):R  
 c$$$     $     ,PTCLOUD_XZ(NCLOUDS_XZ):I'  
 c$$$c$$$     $     ,'//block6  
 c$$$     $     )  
 c$$$  
 c$$$        
 c$$$      return  
 c$$$      end  
 ***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***  
 *  
 *  
 *  
 *  
 *  
 *  
 *  
 *  
 *  
 ***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***  
 c$$$      subroutine book_level2  
 c$$$c*****************************************************  
 c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$cccccc 07/10/2005 modified by elena vannuccini --> (2)  
 c$$$c*****************************************************  
 c$$$  
 c$$$      include 'commontracker.f'  
 c$$$      include 'common_momanhough.f'  
 c$$$      include 'level2.f'  
 c$$$  
 c$$$      character*35 block1,block2  
 c$$$  
 c$$$c      print*,'__________ booking LEVEL2 n-tuple __________'  
 c$$$  
 c$$$*     LEVEL1 ntuple:  
 c$$$      call HBNT(ntp_level2,'LEVEL2',' ')  
 c$$$        
 c$$$c*****************************************************  
 c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$c      call HBNAME(ntp_level2,'EVENT',good2,'GOOD2:L,NEV2:I')  
 c$$$cccccc 06/10/2005 modified by elena vannuccini  
 c$$$c      call HBNAME(ntp_level2,'GENERAL',good2,'GOOD2:L,NEV2:I  
 c$$$c     $     ,WHIC_CALIB:I,SWCODE:I')  
 c$$$      call HBNAME(ntp_level2,'GENERAL',good2,'GOOD2:L,NEV2:I  
 c$$$     $     ,WHICH_CALIB:I,SWCODE:I,CRC(12):L')  
 c$$$c*********************************************************  
 c$$$    
 c$$$  
 c$$$c$$$# ifndef TEST2003  
 c$$$c$$$  
 c$$$c$$$      call HBNAME(ntp_level2,'CPU',pkt_type  
 c$$$c$$$     $     ,'PKT_TYPE:I::[0,50]  
 c$$$c$$$     $     ,PKT_NUM:I  
 c$$$c$$$     $     ,OBT:I'//  
 c$$$c$$$c********************************************************  
 c$$$c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$c$$$c     $     ,WHICH_CALIB:I::[0,50]')  
 c$$$c$$$     $     ',CPU_CRC:L')  
 c$$$c$$$c********************************************************  
 c$$$c$$$  
 c$$$c$$$# endif  
 c$$$  
 c$$$ 417  format('NTRK:I::[0,',I4,']')  
 c$$$ 418  format(',IMAGE(NTRK):I::[0,',I4,']')  
 c$$$      write(block1,417)NTRKMAX  
 c$$$      write(block2,418)NTRKMAX  
 c$$$      call HBNAME(ntp_level2,'TRACKS',NTRK,  
 c$$$     $     block1//  
 c$$$     $     block2//'  
 c$$$     $     ,XM(6,NTRK):R  
 c$$$     $     ,YM(6,NTRK):R  
 c$$$     $     ,ZM(6,NTRK):R  
 c$$$     $     ,RESX(6,NTRK):R  
 c$$$     $     ,RESY(6,NTRK):R  
 c$$$     $     ,AL(5,NTRK):R  
 c$$$     $     ,COVAL(5,5,NTRK):R  
 c$$$     $     ,CHI2(NTRK):R  
 c$$$     $     ,XGOOD(6,NTRK):I::[0,1]  
 c$$$     $     ,YGOOD(6,NTRK):I::[0,1]  
 c$$$     $     ,XV(6,NTRK):R  
 c$$$     $     ,YV(6,NTRK):R  
 c$$$     $     ,ZV(6,NTRK):R  
 c$$$     $     ,AXV(6,NTRK):R  
 c$$$     $     ,AYV(6,NTRK):R'//  
 c$$$c*****************************************************  
 c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$c     $     ,DEDXP(6,NTRK):R'//  
 c$$$c     $     ')  
 c$$$     $     ',DEDX_X(6,NTRK):R  
 c$$$     $     ,DEDX_Y(6,NTRK):R'//  
 c$$$c****************************************************  
 c$$$cccccc 06/10/2005 modified by elena vannuccini  
 c$$$c     $     ,CRC(12):L  
 c$$$     $     ',BdL(NTRK):R'  
 c$$$     $     )  
 c$$$c****************************************************  
 c$$$  
 c$$$    
 c$$$      call HBNAME(ntp_level2,'SINGLETX',nclsx,  
 c$$$c*****************************************************  
 c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$c     $     'NCLSX(6):I,NCLSY(6):I')  
 c$$$     $     'NCLSX:I::[0,500],PLANEX(NCLSX):I  
 c$$$     $     ,XS(2,NCLSX):R,SGNLXS(NCLSX):R') !(2)  
 c$$$c    $     ,XS(NCLSX):R,SGNLXS(NCLSX):R')   !(2)  
 c$$$      call HBNAME(ntp_level2,'SINGLETY',nclsy,  
 c$$$     $     'NCLSY:I::[0,500],PLANEY(NCLSY):I  
 c$$$     $     ,YS(2,NCLSY):R,SGNLYS(NCLSY):R') !(2)  
 c$$$c    $     ,YS(NCLSY):R,SGNLYS(NCLSY):R')   !(2)  
 c$$$      return  
 c$$$      end  
   
 c$$$      subroutine fill_level2_clouds  
 c$$$c*****************************************************  
 c$$$c     29/11/2005 created by elena vannuccini  
 c$$$c*****************************************************  
 c$$$  
 c$$$*     -------------------------------------------------------  
 c$$$*     This routine fills the  variables related to the hough  
 c$$$*     transform, for the debig n-tuple  
 c$$$*     -------------------------------------------------------  
 c$$$  
 c$$$      include 'commontracker.f'  
 c$$$      include 'common_momanhough.f'  
 c$$$      include 'common_level2debug.f'  
 c$$$      include 'level2.f'  
 c$$$  
 c$$$      good2_nt=.true.!good2  
 c$$$c      nev2_nt=nev2  
 c$$$        
 c$$$      if(.false.  
 c$$$     $     .or.ntrpt.gt.ntrpt_max_nt  
 c$$$     $     .or.ndblt.gt.ndblt_max_nt  
 c$$$     $     .or.NCLOUDS_XZ.gt.ncloxz_max  
 c$$$     $     .or.NCLOUDS_yZ.gt.ncloyz_max  
 c$$$     $     )then  
 c$$$         good2_nt=.false.  
 c$$$         ntrpt_nt=0  
 c$$$         ndblt_nt=0  
 c$$$         NCLOUDS_XZ_nt=0  
 c$$$         NCLOUDS_YZ_nt=0  
 c$$$      else  
 c$$$         ndblt_nt=ndblt  
 c$$$         ntrpt_nt=ntrpt  
 c$$$         if(ndblt.ne.0)then  
 c$$$            do id=1,ndblt  
 c$$$               alfayz1_nt(id)=alfayz1(id) !Y0  
 c$$$               alfayz2_nt(id)=alfayz2(id) !tg theta-yz  
 c$$$c               db_cloud_nt(id)=db_cloud(id)  
 c$$$            enddo  
 c$$$         endif  
 c$$$         if(ndblt.ne.0)then  
 c$$$            do it=1,ntrpt  
 c$$$               alfaxz1_nt(it)=alfaxz1(it) !X0  
 c$$$               alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz  
 c$$$               alfaxz3_nt(it)=alfaxz3(it) !1/r  
 c$$$c               tr_cloud_nt(it)=tr_cloud(it)  
 c$$$            enddo  
 c$$$         endif  
 c$$$         nclouds_yz_nt=nclouds_yz  
 c$$$         nclouds_xz_nt=nclouds_xz  
 c$$$         if(nclouds_yz.ne.0)then  
 c$$$            nnn=0  
 c$$$            do iyz=1,nclouds_yz  
 c$$$               ptcloud_yz_nt(iyz)=ptcloud_yz(iyz)  
 c$$$               alfayz1_av_nt(iyz)=alfayz1_av(iyz)  
 c$$$               alfayz2_av_nt(iyz)=alfayz2_av(iyz)  
 c$$$               nnn=nnn+ptcloud_yz(iyz)  
 c$$$            enddo  
 c$$$            do ipt=1,nnn  
 c$$$               db_cloud_nt(ipt)=db_cloud(ipt)  
 c$$$            enddo  
 c$$$c            print*,'#### ntupla #### ptcloud_yz '  
 c$$$c     $           ,(ptcloud_yz(i),i=1,nclouds_yz)  
 c$$$c            print*,'#### ntupla #### db_cloud '  
 c$$$c     $           ,(db_cloud(i),i=1,nnn)  
 c$$$         endif  
 c$$$         if(nclouds_xz.ne.0)then  
 c$$$            nnn=0  
 c$$$            do ixz=1,nclouds_xz  
 c$$$               ptcloud_xz_nt(ixz)=ptcloud_xz(ixz)  
 c$$$               alfaxz1_av_nt(ixz)=alfaxz1_av(ixz)  
 c$$$               alfaxz2_av_nt(ixz)=alfaxz2_av(ixz)  
 c$$$               alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)  
 c$$$               nnn=nnn+ptcloud_xz(ixz)                
 c$$$            enddo  
 c$$$            do ipt=1,nnn  
 c$$$               tr_cloud_nt(ipt)=tr_cloud(ipt)  
 c$$$            enddo  
 c$$$c            print*,'#### ntupla #### ptcloud_xz '  
 c$$$c     $           ,(ptcloud_xz(i),i=1,nclouds_xz)  
 c$$$c            print*,'#### ntupla #### tr_cloud '  
 c$$$c     $           ,(tr_cloud(i),i=1,nnn)  
 c$$$         endif  
 c$$$      endif  
 c$$$      end  
1339                
1340    
1341  ***************************************************  ***************************************************
# Line 1769  c      include 'level1.f' Line 1393  c      include 'level1.f'
1393        enddo        enddo
1394  *     mask views with too many clusters  *     mask views with too many clusters
1395        do iv=1,nviews        do iv=1,nviews
1396           if( ncl_view(iv).gt. nclustermax)then           if( ncl_view(iv).gt. nclusterlimit)then
1397              mask_view(iv) = 1              mask_view(iv) = 1
1398              if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '              if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '
1399       $           ,nclustermax,' on view ', iv,' --> masked!'       $           ,nclusterlimit,' on view ', iv,' --> masked!'
1400           endif           endif
1401        enddo        enddo
1402    
# Line 1789  c      include 'level1.f' Line 1413  c      include 'level1.f'
1413  *     ----------------------------------------------------  *     ----------------------------------------------------
1414  *     cut on charge (X VIEW)  *     cut on charge (X VIEW)
1415  *     ----------------------------------------------------  *     ----------------------------------------------------
1416           if(dedx(icx).lt.dedx_x_min)then           if(sgnl(icx).lt.dedx_x_min)then
1417              cl_single(icx)=0              cl_single(icx)=0
1418              goto 10              goto 10
1419           endif           endif
# Line 1839  c     endif Line 1463  c     endif
1463  *     ----------------------------------------------------  *     ----------------------------------------------------
1464  *     cut on charge (Y VIEW)  *     cut on charge (Y VIEW)
1465  *     ----------------------------------------------------  *     ----------------------------------------------------
1466              if(dedx(icy).lt.dedx_y_min)then              if(sgnl(icy).lt.dedx_y_min)then
1467                 cl_single(icy)=0                 cl_single(icy)=0
1468                 goto 20                 goto 20
1469              endif              endif
# Line 1885  c     endif Line 1509  c     endif
1509  *     charge correlation  *     charge correlation
1510  *     (modified to be applied only below saturation... obviously)  *     (modified to be applied only below saturation... obviously)
1511    
1512                 if(  .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx)                 if(  .not.(sgnl(icy).gt.chsaty.and.sgnl(icx).gt.chsatx)
1513       $              .and.       $              .and.
1514       $              .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx)       $              .not.(sgnl(icy).lt.chmipy.and.sgnl(icx).lt.chmipx)
1515       $              .and.       $              .and.
1516       $              (badclx.eq.1.and.badcly.eq.1)       $              (badclx.eq.1.and.badcly.eq.1)
1517       $              .and.       $              .and.
1518       $              .true.)then       $              .true.)then
1519    
1520                    ddd=(dedx(icy)                    ddd=(sgnl(icy)
1521       $                 -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))       $                 -kch(nplx,nldx)*sgnl(icx)-cch(nplx,nldx))
1522                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
1523    
1524  c                  cut = chcut * sch(nplx,nldx)  c                  cut = chcut * sch(nplx,nldx)
1525    
1526                    sss=(kch(nplx,nldx)*dedx(icy)+dedx(icx)                    sss=(kch(nplx,nldx)*sgnl(icy)+sgnl(icx)
1527       $                 -kch(nplx,nldx)*cch(nplx,nldx))       $                 -kch(nplx,nldx)*cch(nplx,nldx))
1528                    sss=sss/sqrt(kch(nplx,nldx)**2+1)                    sss=sss/sqrt(kch(nplx,nldx)**2+1)
1529                    cut = chcut * (16 + sss/50.)                    cut = chcut * (16 + sss/50.)
# Line 1908  c                  cut = chcut * sch(npl Line 1532  c                  cut = chcut * sch(npl
1532                       goto 20    !charge not consistent                       goto 20    !charge not consistent
1533                    endif                    endif
1534                 endif                 endif
1535                  
1536  *     ------------------> COUPLE <------------------                 if(ncp_plane(nplx).gt.ncouplemax)then
 *     check to do not overflow vector dimentions  
 c$$$               if(ncp_plane(nplx).gt.ncouplemax)then  
 c$$$                  if(DEBUG)print*,  
 c$$$     $                    ' ** warning ** number of identified'//  
 c$$$     $                    ' couples on plane ',nplx,  
 c$$$     $                    ' exceeds vector dimention'//  
 c$$$     $                    ' ( ',ncouplemax,' )'  
 c$$$c     good2=.false.  
 c$$$c     goto 880   !fill ntp and go to next event  
 c$$$                  iflag=1  
 c$$$                  return  
 c$$$               endif  
                 
                if(ncp_plane(nplx).eq.ncouplemax)then  
1537                    if(verbose)print*,                    if(verbose)print*,
1538       $                 '** warning ** number of identified '//       $                 '** warning ** number of identified '//
1539       $                 'couples on plane ',nplx,       $                 'couples on plane ',nplx,
1540       $                 'exceeds vector dimention '       $                 'exceeds vector dimention '
1541       $                 ,'( ',ncouplemax,' )'       $                 ,'( ',ncouplemax,' ) --> masked!'
 c     good2=.false.  
 c     goto 880   !fill ntp and go to next event  
1542                    mask_view(nviewx(nplx)) = 2                    mask_view(nviewx(nplx)) = 2
1543                    mask_view(nviewy(nply)) = 2                    mask_view(nviewy(nply)) = 2
1544  c                  iflag=1                    goto 10
 c                  return  
1545                 endif                 endif
1546                                
1547    *     ------------------> COUPLE <------------------
1548                 ncp_plane(nplx) = ncp_plane(nplx) + 1                 ncp_plane(nplx) = ncp_plane(nplx) + 1
1549                 clx(nplx,ncp_plane(nplx))=icx                 clx(nplx,ncp_plane(nplx))=icx
1550                 cly(nply,ncp_plane(nplx))=icy                 cly(nply,ncp_plane(nplx))=icy
1551                 cl_single(icx)=0                 cl_single(icx)=0
1552                 cl_single(icy)=0                 cl_single(icy)=0
             endif                                
1553  *     ----------------------------------------------  *     ----------------------------------------------
1554    
1555                endif                              
1556    
1557   20         continue   20         continue
1558           enddo                  !end loop on clusters(Y)           enddo                  !end loop on clusters(Y)
1559                    
# Line 1971  c                  return Line 1580  c                  return
1580        do ip=1,6        do ip=1,6
1581           ncp_tot = ncp_tot + ncp_plane(ip)           ncp_tot = ncp_tot + ncp_plane(ip)
1582        enddo        enddo
 c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  
         
 c$$$      if(ncp_tot.gt.ncp_max)then  
 c$$$         if(verbose)print*,  
 c$$$     $           '** warning ** number of identified '//  
 c$$$     $           'couples exceeds upper limit for Hough tr. '  
 c$$$     $           ,'( ',ncp_max,' )'              
 c$$$         iflag=1  
 c$$$         return  
 c$$$      endif  
1583                
1584        return        return
1585        end        end
# Line 1993  c$$$      endif Line 1592  c$$$      endif
1592  *                                                 *  *                                                 *
1593  *                                                 *  *                                                 *
1594  **************************************************  **************************************************
 c$$$      subroutine cl_to_couples_nocharge(iflag)  
 c$$$  
 c$$$      include 'commontracker.f'  
 c$$$      include 'level1.f'  
 c$$$      include 'common_momanhough.f'  
 c$$$c      include 'momanhough_init.f'  
 c$$$      include 'calib.f'  
 c$$$c      include 'level1.f'  
 c$$$  
 c$$$  
 c$$$*     output flag  
 c$$$*     --------------  
 c$$$*     0 = good event  
 c$$$*     1 = bad event  
 c$$$*     --------------  
 c$$$      integer iflag  
 c$$$  
 c$$$      integer badseed,badcl  
 c$$$  
 c$$$*     init variables  
 c$$$      ncp_tot=0  
 c$$$      do ip=1,nplanes  
 c$$$         do ico=1,ncouplemax  
 c$$$            clx(ip,ico)=0  
 c$$$            cly(ip,ico)=0  
 c$$$         enddo  
 c$$$         ncp_plane(ip)=0  
 c$$$         do icl=1,nclstrmax_level2  
 c$$$            cls(ip,icl)=1  
 c$$$         enddo  
 c$$$         ncls(ip)=0  
 c$$$      enddo  
 c$$$      do icl=1,nclstrmax_level2  
 c$$$         cl_single(icl)=1  
 c$$$         cl_good(icl)=0  
 c$$$      enddo  
 c$$$        
 c$$$*     start association  
 c$$$      ncouples=0  
 c$$$      do icx=1,nclstr1          !loop on cluster (X)  
 c$$$         if(mod(VIEW(icx),2).eq.1)goto 10  
 c$$$          
 c$$$*     ----------------------------------------------------  
 c$$$*     cut on charge (X VIEW)  
 c$$$         if(dedx(icx).lt.dedx_x_min)then  
 c$$$            cl_single(icx)=0  
 c$$$            goto 10  
 c$$$         endif  
 c$$$*     cut BAD (X VIEW)              
 c$$$         badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))  
 c$$$         ifirst=INDSTART(icx)  
 c$$$         if(icx.ne.nclstr1) then  
 c$$$            ilast=INDSTART(icx+1)-1  
 c$$$         else  
 c$$$            ilast=TOTCLLENGTH  
 c$$$         endif  
 c$$$         badcl=badseed  
 c$$$         do igood=-ngoodstr,ngoodstr  
 c$$$            ibad=1  
 c$$$            if((INDMAX(icx)+igood).gt.ifirst.and.  
 c$$$     $           (INDMAX(icx)+igood).lt.ilast.and.  
 c$$$     $           .true.)then  
 c$$$               ibad=BAD(VIEW(icx),  
 c$$$     $              nvk(MAXS(icx)+igood),  
 c$$$     $              nst(MAXS(icx)+igood))  
 c$$$            endif  
 c$$$            badcl=badcl*ibad  
 c$$$         enddo  
 c$$$         if(badcl.eq.0)then     !<<<<<<<<<<<<<< BAD cut  
 c$$$            cl_single(icx)=0    !<<<<<<<<<<<<<< BAD cut  
 c$$$            goto 10             !<<<<<<<<<<<<<< BAD cut  
 c$$$         endif                  !<<<<<<<<<<<<<< BAD cut  
 c$$$*     ----------------------------------------------------  
 c$$$          
 c$$$         cl_good(icx)=1  
 c$$$         nplx=npl(VIEW(icx))  
 c$$$         nldx=nld(MAXS(icx),VIEW(icx))  
 c$$$          
 c$$$         do icy=1,nclstr1       !loop on cluster (Y)  
 c$$$            if(mod(VIEW(icy),2).eq.0)goto 20  
 c$$$              
 c$$$*     ----------------------------------------------------  
 c$$$*     cut on charge (Y VIEW)  
 c$$$            if(dedx(icy).lt.dedx_y_min)then  
 c$$$               cl_single(icy)=0  
 c$$$               goto 20  
 c$$$            endif  
 c$$$*     cut BAD (Y VIEW)              
 c$$$            badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))  
 c$$$            ifirst=INDSTART(icy)  
 c$$$            if(icy.ne.nclstr1) then  
 c$$$               ilast=INDSTART(icy+1)-1  
 c$$$            else  
 c$$$               ilast=TOTCLLENGTH  
 c$$$            endif  
 c$$$            badcl=badseed  
 c$$$            do igood=-ngoodstr,ngoodstr  
 c$$$               ibad=1  
 c$$$               if((INDMAX(icy)+igood).gt.ifirst.and.  
 c$$$     $              (INDMAX(icy)+igood).lt.ilast.and.  
 c$$$     $              .true.)  
 c$$$     $              ibad=BAD(VIEW(icy),  
 c$$$     $              nvk(MAXS(icy)+igood),  
 c$$$     $              nst(MAXS(icy)+igood))  
 c$$$               badcl=badcl*ibad  
 c$$$            enddo  
 c$$$            if(badcl.eq.0)then  !<<<<<<<<<<<<<< BAD cut  
 c$$$               cl_single(icy)=0 !<<<<<<<<<<<<<< BAD cut  
 c$$$               goto 20          !<<<<<<<<<<<<<< BAD cut  
 c$$$            endif               !<<<<<<<<<<<<<< BAD cut  
 c$$$*     ----------------------------------------------------  
 c$$$              
 c$$$              
 c$$$            cl_good(icy)=1                    
 c$$$            nply=npl(VIEW(icy))  
 c$$$            nldy=nld(MAXS(icy),VIEW(icy))  
 c$$$              
 c$$$*     ----------------------------------------------  
 c$$$*     CONDITION TO FORM A COUPLE  
 c$$$*     ----------------------------------------------  
 c$$$*     geometrical consistency (same plane and ladder)  
 c$$$            if(nply.eq.nplx.and.nldy.eq.nldx)then  
 c$$$*     charge correlation  
 c$$$*     ===========================================================  
 c$$$*     this version of the subroutine is used for the calibration  
 c$$$*     thus charge-correlation selection is obviously removed  
 c$$$*     ===========================================================  
 c$$$c$$$               ddd=(dedx(icy)  
 c$$$c$$$     $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))  
 c$$$c$$$               ddd=ddd/sqrt(kch(nplx,nldx)**2+1)  
 c$$$c$$$               cut=chcut*sch(nplx,nldx)  
 c$$$c$$$               if(abs(ddd).gt.cut)goto 20 !charge not consistent  
 c$$$*     ===========================================================  
 c$$$                
 c$$$                
 c$$$*     ------------------> COUPLE <------------------  
 c$$$*     check to do not overflow vector dimentions  
 c$$$c$$$               if(ncp_plane(nplx).gt.ncouplemax)then  
 c$$$c$$$                  if(DEBUG)print*,  
 c$$$c$$$     $                    ' ** warning ** number of identified'//  
 c$$$c$$$     $                    ' couples on plane ',nplx,  
 c$$$c$$$     $                    ' exceeds vector dimention'//  
 c$$$c$$$     $                    ' ( ',ncouplemax,' )'  
 c$$$c$$$c     good2=.false.  
 c$$$c$$$c     goto 880   !fill ntp and go to next event  
 c$$$c$$$                  iflag=1  
 c$$$c$$$                  return  
 c$$$c$$$               endif  
 c$$$                
 c$$$               if(ncp_plane(nplx).eq.ncouplemax)then  
 c$$$                  if(verbose)print*,  
 c$$$     $                 '** warning ** number of identified '//  
 c$$$     $                 'couples on plane ',nplx,  
 c$$$     $                 'exceeds vector dimention '  
 c$$$     $                 ,'( ',ncouplemax,' )'  
 c$$$c     good2=.false.  
 c$$$c     goto 880   !fill ntp and go to next event                      
 c$$$                  iflag=1  
 c$$$                  return  
 c$$$               endif  
 c$$$                
 c$$$               ncp_plane(nplx) = ncp_plane(nplx) + 1  
 c$$$               clx(nplx,ncp_plane(nplx))=icx  
 c$$$               cly(nply,ncp_plane(nplx))=icy  
 c$$$               cl_single(icx)=0  
 c$$$               cl_single(icy)=0  
 c$$$            endif                                
 c$$$*     ----------------------------------------------  
 c$$$  
 c$$$ 20         continue  
 c$$$         enddo                  !end loop on clusters(Y)  
 c$$$          
 c$$$ 10      continue  
 c$$$      enddo                     !end loop on clusters(X)  
 c$$$        
 c$$$        
 c$$$      do icl=1,nclstr1  
 c$$$         if(cl_single(icl).eq.1)then  
 c$$$            ip=npl(VIEW(icl))  
 c$$$            ncls(ip)=ncls(ip)+1  
 c$$$            cls(ip,ncls(ip))=icl  
 c$$$         endif  
 c$$$      enddo  
 c$$$        
 c$$$        
 c$$$      if(DEBUG)then  
 c$$$         print*,'clusters  ',nclstr1  
 c$$$         print*,'good    ',(cl_good(i),i=1,nclstr1)  
 c$$$         print*,'singles ',(cl_single(i),i=1,nclstr1)  
 c$$$         print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)  
 c$$$      endif  
 c$$$        
 c$$$      do ip=1,6  
 c$$$         ncp_tot=ncp_tot+ncp_plane(ip)  
 c$$$      enddo  
 c$$$c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  
 c$$$        
 c$$$      if(ncp_tot.gt.ncp_max)then  
 c$$$         if(verbose)print*,  
 c$$$     $           '** warning ** number of identified '//  
 c$$$     $           'couples exceeds upper limit for Hough tr. '  
 c$$$     $           ,'( ',ncp_max,' )'              
 c$$$c            good2=.false.  
 c$$$c     goto 880       !fill ntp and go to next event  
 c$$$         iflag=1  
 c$$$         return  
 c$$$      endif  
 c$$$        
 c$$$      return  
 c$$$      end  
 c$$$  
         
 ***************************************************  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 **************************************************  
1595    
1596        subroutine cp_to_doubtrip(iflag)        subroutine cp_to_doubtrip(iflag)
 c*****************************************************  
 c     02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
1597    
1598        include 'commontracker.f'        include 'commontracker.f'
1599        include 'level1.f'        include 'level1.f'
1600        include 'common_momanhough.f'        include 'common_momanhough.f'
 c      include 'momanhough_init.f'  
1601        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
1602        include 'common_mini_2.f'        include 'common_mini_2.f'
1603        include 'calib.f'        include 'calib.f'
 c      include 'level1.f'  
1604    
1605    
1606  *     output flag  *     output flag
# Line 2260  c      double precision xm3,ym3,zm3 Line 1634  c      double precision xm3,ym3,zm3
1634  *     -----------------------------  *     -----------------------------
1635    
1636    
1637    *     --------------------------------------------
1638    *     put a limit to the maximum number of couples
1639    *     per plane, in order to apply hough transform
1640    *     (couples recovered during track refinement)
1641    *     --------------------------------------------
1642          do ip=1,nplanes
1643             if(ncp_plane(ip).gt.ncouplelimit)then
1644                mask_view(nviewx(ip)) = 8
1645                mask_view(nviewy(ip)) = 8
1646             endif
1647          enddo
1648    
1649    
1650        ndblt=0                   !number of doublets        ndblt=0                   !number of doublets
1651        ntrpt=0                   !number of triplets        ntrpt=0                   !number of triplets
1652                
1653        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1
1654           do is1=1,2             !loop on sensors - COPPIA 1           if(  mask_view(nviewx(ip1)).ne.0 .or.
1655                     $        mask_view(nviewy(ip1)).ne.0 )goto 10 !skip plane
1656             do is1=1,2             !loop on sensors - COPPIA 1            
1657              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1
1658                 icx1=clx(ip1,icp1)                 icx1=clx(ip1,icp1)
1659                 icy1=cly(ip1,icp1)                 icy1=cly(ip1,icp1)
1660  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)
1661                 call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)  c               call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)
1662                   call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.)
1663                 xm1=xPAM                 xm1=xPAM
1664                 ym1=yPAM                 ym1=yPAM
1665                 zm1=zPAM                                   zm1=zPAM                  
1666  c     print*,'***',is1,xm1,ym1,zm1  c     print*,'***',is1,xm1,ym1,zm1
1667    
1668                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
1669                      if(  mask_view(nviewx(ip2)).ne.0 .or.
1670         $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane
1671                    do is2=1,2    !loop on sensors -ndblt COPPIA 2                    do is2=1,2    !loop on sensors -ndblt COPPIA 2
1672                                            
1673                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
# Line 2284  c     print*,'***',is1,xm1,ym1,zm1 Line 1675  c     print*,'***',is1,xm1,ym1,zm1
1675                          icy2=cly(ip2,icp2)                          icy2=cly(ip2,icp2)
1676  c                        call xyz_PAM  c                        call xyz_PAM
1677  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)
1678    c                        call xyz_PAM
1679    c     $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)
1680                          call xyz_PAM                          call xyz_PAM
1681       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.)
1682                          xm2=xPAM                          xm2=xPAM
1683                          ym2=yPAM                          ym2=yPAM
1684                          zm2=zPAM                          zm2=zPAM
# Line 2337  c$$$ Line 1730  c$$$
1730  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1731    
1732    
1733                          if(ip2.eq.nplanes)goto 30 !no possible combination with 3 couples                          if(ip2.eq.nplanes)goto 31 !no possible combination with 3 couples
1734    
1735                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3
1736                               if(  mask_view(nviewx(ip3)).ne.0 .or.
1737         $                          mask_view(nviewy(ip3)).ne.0 )goto 30 !skip plane
1738                             do is3=1,2 !loop on sensors - COPPIA 3                             do is3=1,2 !loop on sensors - COPPIA 3
1739                                                                
1740                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3
# Line 2346  c$$$ Line 1742  c$$$
1742                                   icy3=cly(ip3,icp3)                                   icy3=cly(ip3,icp3)
1743  c                                 call xyz_PAM  c                                 call xyz_PAM
1744  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)
1745    c                                 call xyz_PAM
1746    c     $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)
1747                                   call xyz_PAM                                   call xyz_PAM
1748       $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)       $                                (icx3,icy3,is3,PFAdef,PFAdef
1749         $                                ,0.,0.,0.,0.)
1750                                   xm3=xPAM                                   xm3=xPAM
1751                                   ym3=yPAM                                   ym3=yPAM
1752                                   zm3=zPAM                                   zm3=zPAM
# Line 2415  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp Line 1814  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp
1814                                endif                                endif
1815                             enddo !end loop on COPPIA 3                             enddo !end loop on COPPIA 3
1816                          enddo   !end loop on sensors - COPPIA 3                          enddo   !end loop on sensors - COPPIA 3
1817     30                     continue
1818                       enddo      !end loop on planes  - COPPIA 3                       enddo      !end loop on planes  - COPPIA 3
1819   30                  continue   31                  continue
1820                                            
1821   1                enddo         !end loop on COPPIA 2   1                enddo         !end loop on COPPIA 2
1822                 enddo            !end loop on sensors - COPPIA 2                 enddo            !end loop on sensors - COPPIA 2
1823     20            continue
1824              enddo               !end loop on planes  - COPPIA 2              enddo               !end loop on planes  - COPPIA 2
1825                            
1826           enddo                  !end loop on COPPIA1           enddo                  !end loop on COPPIA1
1827        enddo                     !end loop on sensors - COPPIA 1        enddo                     !end loop on sensors - COPPIA 1
1828     10   continue
1829        enddo                     !end loop on planes  - COPPIA 1        enddo                     !end loop on planes  - COPPIA 1
1830                
1831        if(DEBUG)then        if(DEBUG)then
# Line 2993  c      real fitz(nplanes)        !z coor Line 2395  c      real fitz(nplanes)        !z coor
2395                 nplused=nplused+ hit_plane(ip)                 nplused=nplused+ hit_plane(ip)
2396              enddo              enddo
2397                            
2398              if(nplused.lt.nplxz_min)goto 888 !next doublet  c            if(nplused.lt.nplxz_min)goto 888 !next doublet
2399                if(nplused.lt.nplyz_min)goto 888 !next doublet
2400              if(ncp_ok.lt.ncpok)goto 888 !next cloud              if(ncp_ok.lt.ncpok)goto 888 !next cloud
2401                            
2402              if(DEBUG)then              if(DEBUG)then
# Line 3097  c$$$            if(AL_INI(5).gt.defmax)g Line 2500  c$$$            if(AL_INI(5).gt.defmax)g
2500  *                                   *************************  *                                   *************************
2501  c                                    call xyz_PAM(icx,icy,is,  c                                    call xyz_PAM(icx,icy,is,
2502  c     $                                   'COG2','COG2',0.,0.)  c     $                                   'COG2','COG2',0.,0.)
2503    c                                    call xyz_PAM(icx,icy,is, !(1)
2504    c     $                                   PFAdef,PFAdef,0.,0.) !(1)
2505                                      call xyz_PAM(icx,icy,is, !(1)                                      call xyz_PAM(icx,icy,is, !(1)
2506       $                                   PFAdef,PFAdef,0.,0.) !(1)       $                                   PFAdef,PFAdef,0.,0.,0.,0.)
2507  *                                   *************************  *                                   *************************
2508  *                                   -----------------------------  *                                   -----------------------------
2509                                      xgood(nplanes-ip+1)=1.                                      xgood(nplanes-ip+1)=1.
# Line 3126  c$$$                              enddo Line 2531  c$$$                              enddo
2531                                jstep=0 !number of  minimization steps                                jstep=0 !number of  minimization steps
2532                                iprint=0                                iprint=0
2533  c                              if(DEBUG)iprint=1  c                              if(DEBUG)iprint=1
2534                                if(DEBUG)iprint=1                                if(DEBUG)iprint=2
2535                                call mini2(jstep,ifail,iprint)                                call mini2(jstep,ifail,iprint)
2536                                if(ifail.ne.0) then                                if(ifail.ne.0) then
2537                                   if(DEBUG)then                                   if(DEBUG)then
# Line 3270  c      include 'momanhough_init.f' Line 2675  c      include 'momanhough_init.f'
2675  c      include 'level1.f'  c      include 'level1.f'
2676        include 'calib.f'        include 'calib.f'
2677    
   
2678  *     flag to chose PFA  *     flag to chose PFA
2679        character*10 PFA        character*10 PFA
2680        common/FINALPFA/PFA        common/FINALPFA/PFA
2681    
2682          real xp,yp,zp
2683          real xyzp(3),bxyz(3)
2684          equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3))
2685    
2686  *     =================================================  *     =================================================
2687  *     new estimate of positions using ETA algorithm  *     new estimate of positions using ETA algorithm
2688  *                          and  *                          and
# Line 3283  c      include 'level1.f' Line 2691  c      include 'level1.f'
2691        call track_init        call track_init
2692        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
2693    
2694             xP=XV_STORE(nplanes-ip+1,ibest)
2695             yP=YV_STORE(nplanes-ip+1,ibest)
2696             zP=ZV_STORE(nplanes-ip+1,ibest)
2697             call gufld(xyzp,bxyz)
2698    c$$$         bxyz(1)=0
2699    c$$$         bxyz(2)=0
2700    c$$$         bxyz(3)=0
2701  *     |||||||||||||||||||||||||||||||||||||||||||||||||  *     |||||||||||||||||||||||||||||||||||||||||||||||||
2702  *     -------------------------------------------------  *     -------------------------------------------------
2703  *     If the plane has been already included, it just  *     If the plane has been already included, it just
# Line 3303  c      include 'level1.f' Line 2718  c      include 'level1.f'
2718       $           ,ip_cp(id),ip       $           ,ip_cp(id),ip
2719              icx=clx(ip,icp)              icx=clx(ip,icp)
2720              icy=cly(ip,icp)              icy=cly(ip,icp)
2721    c            call xyz_PAM(icx,icy,is,
2722    c     $           PFA,PFA,
2723    c     $           AXV_STORE(nplanes-ip+1,ibest),
2724    c     $           AYV_STORE(nplanes-ip+1,ibest))
2725              call xyz_PAM(icx,icy,is,              call xyz_PAM(icx,icy,is,
 c     $           'ETA2','ETA2',  
2726       $           PFA,PFA,       $           PFA,PFA,
2727       $           AXV_STORE(nplanes-ip+1,ibest),       $           AXV_STORE(nplanes-ip+1,ibest),
2728       $           AYV_STORE(nplanes-ip+1,ibest))       $           AYV_STORE(nplanes-ip+1,ibest),
2729         $           bxyz(1),
2730         $           bxyz(2)
2731         $           )
2732  c$$$  call xyz_PAM(icx,icy,is,  c$$$  call xyz_PAM(icx,icy,is,
2733  c$$$  $              'COG2','COG2',  c$$$  $              'COG2','COG2',
2734  c$$$  $              0.,  c$$$  $              0.,
# Line 3320  c$$$  $              0.) Line 2741  c$$$  $              0.)
2741              resx(nplanes-ip+1) = resxPAM              resx(nplanes-ip+1) = resxPAM
2742              resy(nplanes-ip+1) = resyPAM              resy(nplanes-ip+1) = resyPAM
2743    
2744  c            dedxtrk(nplanes-ip+1) = (dedx(icx)+dedx(icy))/2. !(1)  c            dedxtrk(nplanes-ip+1) = (sgnl(icx)+sgnl(icy))/2. !(1)
2745              dedxtrk_x(nplanes-ip+1)=dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)              dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
2746              dedxtrk_y(nplanes-ip+1)=dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)              dedxtrk_y(nplanes-ip+1)=sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
2747                            
2748  *     |||||||||||||||||||||||||||||||||||||||||||||||||  *     |||||||||||||||||||||||||||||||||||||||||||||||||
2749  *     -------------------------------------------------  *     -------------------------------------------------
# Line 3337  c            dedxtrk(nplanes-ip+1) = (de Line 2758  c            dedxtrk(nplanes-ip+1) = (de
2758                                
2759  *     --------------------------------------------------------------  *     --------------------------------------------------------------
2760  *     determine which ladder and sensor are intersected by the track  *     determine which ladder and sensor are intersected by the track
2761              xP=XV_STORE(nplanes-ip+1,ibest)  c$$$            xP=XV_STORE(nplanes-ip+1,ibest)
2762              yP=YV_STORE(nplanes-ip+1,ibest)  c$$$            yP=YV_STORE(nplanes-ip+1,ibest)
2763              zP=ZV_STORE(nplanes-ip+1,ibest)  c$$$            zP=ZV_STORE(nplanes-ip+1,ibest)
2764              call whichsensor(ip,xP,yP,nldt,ist)              call whichsensor(ip,xP,yP,nldt,ist)
2765  *     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
2766              if(nldt.eq.0.or.ist.eq.0)goto 133              if(nldt.eq.0.or.ist.eq.0)goto 133
# Line 3377  c     $              cl_used(icy).eq.1.o Line 2798  c     $              cl_used(icy).eq.1.o
2798       $              cl_used(icy).ne.0.or. !or the Y cluster is already used !(3)       $              cl_used(icy).ne.0.or. !or the Y cluster is already used !(3)
2799       $              .false.)goto 1188 !then jump to next couple.       $              .false.)goto 1188 !then jump to next couple.
2800  *            *          
2801    c               call xyz_PAM(icx,icy,ist,
2802    c     $              PFA,PFA,
2803    c     $              AXV_STORE(nplanes-ip+1,ibest),
2804    c     $              AYV_STORE(nplanes-ip+1,ibest))
2805                 call xyz_PAM(icx,icy,ist,                 call xyz_PAM(icx,icy,ist,
2806       $              PFA,PFA,       $              PFA,PFA,
 c     $              'ETA2','ETA2',  
2807       $              AXV_STORE(nplanes-ip+1,ibest),       $              AXV_STORE(nplanes-ip+1,ibest),
2808       $              AYV_STORE(nplanes-ip+1,ibest))       $              AYV_STORE(nplanes-ip+1,ibest),
2809         $              bxyz(1),
2810         $              bxyz(2)
2811         $              )
2812                                
2813                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
2814                 distance = distance / RCHI2_STORE(ibest)!<<< MS                 distance = distance / RCHI2_STORE(ibest)!<<< MS
# Line 3396  c     $              'ETA2','ETA2', Line 2823  c     $              'ETA2','ETA2',
2823                    rymm = resyPAM                    rymm = resyPAM
2824                    distmin = distance                    distmin = distance
2825                    idm = id                                      idm = id                  
2826  c                 dedxmm = (dedx(icx)+dedx(icy))/2. !(1)  c                 dedxmm = (sgnl(icx)+sgnl(icy))/2. !(1)
2827                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
2828                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
2829                 endif                 endif
2830   1188          continue   1188          continue
2831              enddo               !end loop on couples on plane icp              enddo               !end loop on couples on plane icp
# Line 3450  c            if(DEBUG)print*,'>>>> try t Line 2877  c            if(DEBUG)print*,'>>>> try t
2877  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
2878                 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)
2879  *                                              !jump to the Y cluster  *                                              !jump to the Y cluster
2880    c               call xyz_PAM(icx,0,ist,
2881    c     $              PFA,PFA,
2882    c     $              AXV_STORE(nplanes-ip+1,ibest),0.)              
2883                 call xyz_PAM(icx,0,ist,                 call xyz_PAM(icx,0,ist,
 c     $              'ETA2','ETA2',  
2884       $              PFA,PFA,       $              PFA,PFA,
2885       $              AXV_STORE(nplanes-ip+1,ibest),0.)                     $              AXV_STORE(nplanes-ip+1,ibest),0.,
2886         $              bxyz(1),
2887         $              bxyz(2)
2888         $              )              
2889                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
2890                 distance = distance / RCHI2_STORE(ibest)!<<< MS                 distance = distance / RCHI2_STORE(ibest)!<<< MS
2891                 if(DEBUG)print*,'( cl-X ',icx                 if(DEBUG)print*,'( cl-X ',icx
# Line 3469  c     $              'ETA2','ETA2', Line 2901  c     $              'ETA2','ETA2',
2901                    rymm = resyPAM                    rymm = resyPAM
2902                    distmin = distance                    distmin = distance
2903                    iclm = icx                    iclm = icx
2904  c                  dedxmm = dedx(icx) !(1)  c                  dedxmm = sgnl(icx) !(1)
2905                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
2906                    dedxmmy = 0.        !(1)                    dedxmmy = 0.        !(1)
2907                 endif                                   endif                  
2908  11881          continue  11881          continue
# Line 3478  c                  dedxmm = dedx(icx) !( Line 2910  c                  dedxmm = dedx(icx) !(
2910  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
2911                 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)
2912  *                                              !jump to the next couple  *                                              !jump to the next couple
2913    c               call xyz_PAM(0,icy,ist,
2914    c     $              PFA,PFA,
2915    c     $              0.,AYV_STORE(nplanes-ip+1,ibest))
2916                 call xyz_PAM(0,icy,ist,                 call xyz_PAM(0,icy,ist,
 c     $              'ETA2','ETA2',  
2917       $              PFA,PFA,       $              PFA,PFA,
2918       $              0.,AYV_STORE(nplanes-ip+1,ibest))       $              0.,AYV_STORE(nplanes-ip+1,ibest),
2919         $              bxyz(1),
2920         $              bxyz(2)
2921         $              )
2922                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
2923                 distance = distance / RCHI2_STORE(ibest)!<<< MS                 distance = distance / RCHI2_STORE(ibest)!<<< MS
2924                 if(DEBUG)print*,'( cl-Y ',icy                 if(DEBUG)print*,'( cl-Y ',icy
# Line 3497  c     $              'ETA2','ETA2', Line 2934  c     $              'ETA2','ETA2',
2934                    rymm = resyPAM                    rymm = resyPAM
2935                    distmin = distance                    distmin = distance
2936                    iclm = icy                    iclm = icy
2937  c                 dedxmm = dedx(icy)  !(1)  c                 dedxmm = sgnl(icy)  !(1)
2938                    dedxmmx = 0.        !(1)                    dedxmmx = 0.        !(1)
2939                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
2940                 endif                                   endif                  
2941  11882          continue  11882          continue
2942              enddo               !end loop on cluster inside couples              enddo               !end loop on cluster inside couples
2943  *----- single clusters -----------------------------------------------      *----- single clusters -----------------------------------------------  
2944    c            print*,'## ncls(',ip,') ',ncls(ip)
2945              do ic=1,ncls(ip)    !loop on single clusters              do ic=1,ncls(ip)    !loop on single clusters
2946                 icl=cls(ip,ic)                 icl=cls(ip,ic)
2947    c              print*,'## ic ',ic,' ist ',ist
2948  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
2949                 if(cl_used(icl).ne.0.or.     !if the cluster is already used !(3)                 if(cl_used(icl).ne.0.or.     !if the cluster is already used !(3)
2950       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match
2951       $              .false.)goto 18882      !jump to the next singlet       $              .false.)goto 18882      !jump to the next singlet
2952                 if(mod(VIEW(icl),2).eq.0)then!<---- X view                 if(mod(VIEW(icl),2).eq.0)then!<---- X view
2953    c                  call xyz_PAM(icl,0,ist,
2954    c     $                 PFA,PFA,
2955    c     $                 AXV_STORE(nplanes-ip+1,ibest),0.)
2956                    call xyz_PAM(icl,0,ist,                    call xyz_PAM(icl,0,ist,
 c     $                 'ETA2','ETA2',  
2957       $                 PFA,PFA,       $                 PFA,PFA,
2958       $                 AXV_STORE(nplanes-ip+1,ibest),0.)       $                 AXV_STORE(nplanes-ip+1,ibest),0.,
2959         $                 bxyz(1),
2960         $                 bxyz(2)
2961         $                 )
2962                 else                         !<---- Y view                 else                         !<---- Y view
2963    c                  call xyz_PAM(0,icl,ist,
2964    c     $                 PFA,PFA,
2965    c     $                 0.,AYV_STORE(nplanes-ip+1,ibest))
2966                    call xyz_PAM(0,icl,ist,                    call xyz_PAM(0,icl,ist,
 c     $                 'ETA2','ETA2',  
2967       $                 PFA,PFA,       $                 PFA,PFA,
2968       $                 0.,AYV_STORE(nplanes-ip+1,ibest))       $                 0.,AYV_STORE(nplanes-ip+1,ibest),
2969         $                 bxyz(1),
2970         $                 bxyz(2)
2971         $                 )
2972                 endif                 endif
2973    
2974                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
2975                 distance = distance / RCHI2_STORE(ibest)!<<< MS                 distance = distance / RCHI2_STORE(ibest)!<<< MS
2976                 if(DEBUG)print*,'( cl-s ',icl                 if(DEBUG)print*,'( cl-s ',icl
2977       $              ,' ) normalized distance ',distance       $              ,' ) normalized distance ',distance,'<',distmin,' ?'
2978                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
2979                      if(DEBUG)print*,'YES'
2980                    xmm_A = xPAM_A                    xmm_A = xPAM_A
2981                    ymm_A = yPAM_A                    ymm_A = yPAM_A
2982                    zmm_A = zPAM_A                    zmm_A = zPAM_A
# Line 3537  c     $                 'ETA2','ETA2', Line 2987  c     $                 'ETA2','ETA2',
2987                    rymm = resyPAM                    rymm = resyPAM
2988                    distmin = distance                      distmin = distance  
2989                    iclm = icl                    iclm = icl
2990  c                  dedxmm = dedx(icl)                   !(1)  c                  dedxmm = sgnl(icl)                   !(1)
2991                    if(mod(VIEW(icl),2).eq.0)then !<---- X view                    if(mod(VIEW(icl),2).eq.0)then !<---- X view
2992                       dedxmmx = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmx = sgnl(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)
2993                       dedxmmy = 0.                       !(1)                       dedxmmy = 0.                       !(1)
2994                    else          !<---- Y view                    else          !<---- Y view
2995                       dedxmmx = 0.                       !(1)                       dedxmmx = 0.                       !(1)
2996                       dedxmmy = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmy = sgnl(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)
2997                    endif                    endif
2998                 endif                                   endif                  
2999  18882          continue  18882          continue
3000              enddo               !end loop on single clusters              enddo               !end loop on single clusters
3001    c            print*,'## distmin ', distmin,' clinc ',clinc
3002              if(distmin.le.clinc)then                                if(distmin.le.clinc)then                  
3003                                
3004                 CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<                     CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<    
# Line 3684  c               endif Line 3134  c               endif
3134    
3135    
3136    
 c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  
 c$$$      real function fbad_cog(ncog,ic)  
 c$$$  
 c$$$  
 c$$$      include 'commontracker.f'  
 c$$$      include 'level1.f'  
 c$$$      include 'calib.f'  
 c$$$  
 c$$$*     --> signal of the central strip  
 c$$$      sc = CLSIGNAL(INDMAX(ic)) !center  
 c$$$  
 c$$$*     signal of adjacent strips  
 c$$$*     --> left  
 c$$$      sl1 = 0                  !left 1  
 c$$$      if(  
 c$$$     $     (INDMAX(ic)-1).ge.INDSTART(ic)  
 c$$$     $     )  
 c$$$     $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))  
 c$$$  
 c$$$      sl2 = 0                  !left 2  
 c$$$      if(  
 c$$$     $     (INDMAX(ic)-2).ge.INDSTART(ic)  
 c$$$     $     )  
 c$$$     $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))  
 c$$$  
 c$$$*     --> right  
 c$$$      sr1 = 0                  !right 1  
 c$$$      if(  
 c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))  
 c$$$     $     .or.  
 c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)  
 c$$$     $     )  
 c$$$     $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))  
 c$$$  
 c$$$      sr2 = 0                  !right 2  
 c$$$      if(  
 c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))  
 c$$$     $     .or.  
 c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)  
 c$$$     $     )  
 c$$$     $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))  
 c$$$  
 c$$$  
 c$$$      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view  
 c$$$         f  = 4.  
 c$$$         si = 8.4  
 c$$$      else                              !X-view  
 c$$$         f  = 6.  
 c$$$         si = 3.9  
 c$$$      endif  
 c$$$  
 c$$$      fbad_cog = 1.  
 c$$$      f0 = 1  
 c$$$      f1 = 1  
 c$$$      f2 = 1  
 c$$$      f3 = 1    
 c$$$      if(sl1.gt.sr1.and.sl1.gt.0.)then  
 c$$$          
 c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f  
 c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f  
 c$$$c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f  
 c$$$  
 c$$$         if(ncog.eq.2.and.sl1.ne.0)then  
 c$$$            fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.)  
 c$$$         elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then  
 c$$$            fbad_cog = 1.  
 c$$$         elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then  
 c$$$            fbad_cog = 1.  
 c$$$         else  
 c$$$            fbad_cog = 1.  
 c$$$         endif  
 c$$$          
 c$$$      elseif(sl1.le.sr1.and.sr1.gt.0.)then  
 c$$$  
 c$$$  
 c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f  
 c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f  
 c$$$c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f  
 c$$$  
 c$$$         if(ncog.eq.2.and.sr1.ne.0)then  
 c$$$            fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.)  
 c$$$         elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then  
 c$$$            fbad_cog = 1.  
 c$$$         elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then  
 c$$$            fbad_cog = 1.  
 c$$$         else  
 c$$$            fbad_cog = 1.  
 c$$$         endif  
 c$$$  
 c$$$      endif  
 c$$$  
 c$$$      fbad_cog = sqrt(fbad_cog)  
 c$$$  
 c$$$      return  
 c$$$      end  
 c$$$  
   
3137    
3138    
3139  *     ****************************************************  *     ****************************************************
# Line 3865  c      include 'level1.f' Line 3218  c      include 'level1.f'
3218           alfayz1_nt(idb)=0                 alfayz1_nt(idb)=0      
3219           alfayz2_nt(idb)=0                 alfayz2_nt(idb)=0      
3220        enddo        enddo
3221        do itr=1,ntrpl_max_nt        do itr=1,ntrpt_max_nt
3222           tr_cloud_nt(itr)=0           tr_cloud_nt(itr)=0
3223           alfaxz1_nt(itr)=0                 alfaxz1_nt(itr)=0      
3224           alfaxz2_nt(itr)=0                 alfaxz2_nt(itr)=0      
# Line 3894  c      include 'level1.f' Line 3247  c      include 'level1.f'
3247          alfayz1(idb)=0                    alfayz1(idb)=0          
3248          alfayz2(idb)=0                    alfayz2(idb)=0          
3249        enddo                            enddo                    
3250        do itr=1,ntrpl_max                do itr=1,ntrpt_max        
3251          tr_cloud(itr)=0                  tr_cloud(itr)=0        
3252          cpxz1(itr)=0                      cpxz1(itr)=0            
3253          cpxz2(itr)=0                      cpxz2(itr)=0            
# Line 3985  c      include 'level1.f' Line 3338  c      include 'level1.f'
3338           yv_nt(ip,ntr)    = sngl(yv(ip))           yv_nt(ip,ntr)    = sngl(yv(ip))
3339           zv_nt(ip,ntr)    = sngl(zv(ip))           zv_nt(ip,ntr)    = sngl(zv(ip))
3340           axv_nt(ip,ntr)   = sngl(axv(ip))           axv_nt(ip,ntr)   = sngl(axv(ip))
3341           ayv_nt(ip,ntr)   = sngl(ayv(ip))           ayv_nt(ip,ntr)   = sngl(ayv(ip))  
3342           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)  c     l'avevo dimenticato!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3343           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)             factor = sqrt(
3344         $        tan( acos(-1.) * sngl(axv(ip)) /180. )**2 +
3345         $        tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 +
3346         $        1. )
3347    c     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3348             dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)
3349             dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)  
3350        
3351           id  = CP_STORE(ip,IDCAND)           id  = CP_STORE(ip,IDCAND)
3352           icl = CLS_STORE(ip,IDCAND)           icl = CLS_STORE(ip,IDCAND)
# Line 4037  c      good2=1!.true. Line 3396  c      good2=1!.true.
3396              if(mod(VIEW(icl),2).eq.0)then !=== X views              if(mod(VIEW(icl),2).eq.0)then !=== X views
3397                 nclsx = nclsx + 1                 nclsx = nclsx + 1
3398                 planex(nclsx) = ip                 planex(nclsx) = ip
3399                 sgnlxs(nclsx) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))!(2)
3400                 clsx(nclsx)   = icl                 clsx(nclsx)   = icl
3401                 do is=1,2                 do is=1,2
3402  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
3403                    call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)  c                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
3404                      call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.,0.,0.)
3405                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2
3406                 enddo                 enddo
3407  c$$$               print*,'nclsx         ',nclsx  c$$$               print*,'nclsx         ',nclsx
# Line 4052  c$$$               print*,'xs(2,nclsx)   Line 3412  c$$$               print*,'xs(2,nclsx)  
3412              else                          !=== Y views              else                          !=== Y views
3413                 nclsy = nclsy + 1                 nclsy = nclsy + 1
3414                 planey(nclsy) = ip                 planey(nclsy) = ip
3415                 sgnlys(nclsy) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))!(2)
3416                 clsy(nclsy)   = icl                 clsy(nclsy)   = icl
3417                 do is=1,2                 do is=1,2
3418  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
3419                    call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)  c                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
3420                      call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.,0.,0.)
3421                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2
3422                 enddo                 enddo
3423  c$$$               print*,'nclsy         ',nclsy  c$$$               print*,'nclsy         ',nclsy
Legend:
Removed from v.1.13  
changed lines
  Added in v.1.19
  ViewVC Help
Powered by ViewVC 1.1.23