/[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.37 by pam-fi, Fri Dec 5 08:26:47 2008 UTC revision 1.44 by pamelats, Wed Dec 30 14:22:23 2009 UTC
# Line 271  c      include 'momanhough_init.f' Line 271  c      include 'momanhough_init.f'
271  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
272  ccc         ntrk=0                 !counter of identified physical tracks  ccc         ntrk=0                 !counter of identified physical tracks
273    
274  11111    continue               !<<<<<<< come here when performing a new search  c11111    continue               !<<<<<<< come here when performing a new search
275             continue                  !<<<<<<< come here when performing a new search
276    
277           if(nclouds_xz.eq.0)goto 880 !go to next event               if(nclouds_xz.eq.0)goto 880 !go to next event    
278           if(nclouds_yz.eq.0)goto 880 !go to next event               if(nclouds_yz.eq.0)goto 880 !go to next event    
# Line 302  c$$$         if(ibest.eq.0)goto 880 !>> Line 303  c$$$         if(ibest.eq.0)goto 880 !>>
303  *     1st) decreasing n.points  *     1st) decreasing n.points
304  *     2nd) increasing chi**2  *     2nd) increasing chi**2
305  *     -------------------------------------------------------  *     -------------------------------------------------------
306           rchi2best=1000000000.           rchi2best=1000000000.
307           ndofbest=0                       ndofbest=0            
308           do i=1,ntracks           do i=1,ntracks
309             ndof=0                           ndof=0              
# Line 376  c         if(DEBUG.EQ.1)iprint=1 Line 377  c         if(DEBUG.EQ.1)iprint=1
377           if(VERBOSE.EQ.1)iprint=1           if(VERBOSE.EQ.1)iprint=1
378           if(DEBUG.EQ.1)iprint=2           if(DEBUG.EQ.1)iprint=2
379           call mini2(jstep,ifail,iprint)           call mini2(jstep,ifail,iprint)
380           if(ifail.ne.0) then  cc         if(ifail.ne.0) then
381             if(ifail.ne.0.or.CHI2.ne.CHI2) then !new
382                if(CHI2.ne.CHI2)CHI2=-9999. !new
383              if(VERBOSE.EQ.1)then              if(VERBOSE.EQ.1)then
384                 print *,                 print *,
385       $              '*** MINIMIZATION FAILURE *** (after refinement) '       $              '*** MINIMIZATION FAILURE *** (after refinement) '
386       $              ,iev       $              ,iev
   
387              endif              endif
388           endif           endif
389                    
# Line 526  c         rchi2=chi2/dble(ndof) Line 528  c         rchi2=chi2/dble(ndof)
528    
529    
530  *     --- and store the results --------------------------------  *     --- and store the results --------------------------------
          ntrk = ntrk + 1                   !counter of found tracks  
          if(.not.FIMAGE  
      $        .and.iimage.eq.0) image(ntrk)= 0  
          if(.not.FIMAGE  
      $        .and.iimage.ne.0)image(ntrk)=ntrk+1 !this is the image of the next  
          if(FIMAGE)     image(ntrk)=ntrk-1 !this is the image of the previous  
          call fill_level2_tracks(ntrk)     !==> good2=.true.  
   
531           if(ntrk.eq.NTRKMAX)then           if(ntrk.eq.NTRKMAX)then
532              if(verbose.eq.1)              if(verbose.eq.1)
533       $           print*,       $           print*,
# Line 543  c         rchi2=chi2/dble(ndof) Line 537  c         rchi2=chi2/dble(ndof)
537  cc            good2=.false.  cc            good2=.false.
538              goto 880            !fill ntp and go to next event              goto 880            !fill ntp and go to next event
539           endif           endif
540    
541             ntrk = ntrk + 1                   !counter of found tracks
542             if(.not.FIMAGE
543         $        .and.iimage.eq.0) image(ntrk)= 0
544             if(.not.FIMAGE
545         $        .and.iimage.ne.0)image(ntrk)=ntrk+1 !this is the image of the next
546             if(FIMAGE)     image(ntrk)=ntrk-1 !this is the image of the previous
547             call fill_level2_tracks(ntrk)     !==> good2=.true.
548    
549    c$$$         if(ntrk.eq.NTRKMAX)then
550    c$$$            if(verbose.eq.1)
551    c$$$     $           print*,
552    c$$$     $           '** warning ** number of identified '//
553    c$$$     $           'tracks exceeds vector dimension '
554    c$$$     $           ,'( ',NTRKMAX,' )'
555    c$$$cc            good2=.false.
556    c$$$            goto 880            !fill ntp and go to next event
557    c$$$         endif
558           if(iimage.ne.0)then           if(iimage.ne.0)then
559              FIMAGE=.true.       !              FIMAGE=.true.       !
560              goto 1212           !>>> fit image-track              goto 1212           !>>> fit image-track
# Line 704  c         resxPAM = RESXAV Line 716  c         resxPAM = RESXAV
716           stripx  = stripx + corr           stripx  = stripx + corr
717           resxPAM = res           resxPAM = res
718    
719   10   endif   10   continue
720          endif
721            
722  *     -----------------  *     -----------------
723  *     CLUSTER Y  *     CLUSTER Y
# Line 750  c         resyPAM = RESYAV Line 763  c         resyPAM = RESYAV
763           stripy  = stripy + corr           stripy  = stripy + corr
764           resyPAM = res           resyPAM = res
765    
766   20   endif   20   continue
767          endif
768    
769    
770  c===========================================================  c===========================================================
# Line 1580  c      include 'level1.f' Line 1594  c      include 'level1.f'
1594        integer iflag        integer iflag
1595    
1596        integer badseed,badclx,badcly        integer badseed,badclx,badcly
1597        
1598          iflag = iflag
1599        if(DEBUG.EQ.1)print*,'cl_to_couples:'        if(DEBUG.EQ.1)print*,'cl_to_couples:'
1600    
1601  cc      if(RECOVER_SINGLETS.and..not.SECOND_SEARCH)goto 80  cc      if(RECOVER_SINGLETS.and..not.SECOND_SEARCH)goto 80
# Line 1766  c                  cut = chcut * sch(npl Line 1781  c                  cut = chcut * sch(npl
1781       $                 ,'( ',ncouplemax,' ) --> masked!'       $                 ,'( ',ncouplemax,' ) --> masked!'
1782  c                  mask_view(nviewx(nplx)) = 2  c                  mask_view(nviewx(nplx)) = 2
1783  c                  mask_view(nviewy(nply)) = 2  c                  mask_view(nviewy(nply)) = 2
1784                    mask_view(nviewx(nplx))= mask_view(nviewx(nplx))+ 2**1                   mask_view(nviewx(nplx))= mask_view(nviewx(nplx))+ 2**1
1785                    mask_view(nviewy(nply))= mask_view(nviewy(nply))+ 2**1                   mask_view(nviewy(nply))= mask_view(nviewy(nply))+ 2**1
1786                    goto 10                    goto 10
1787                 endif                 endif
1788                                
# Line 1795  c                  mask_view(nviewy(nply Line 1810  c                  mask_view(nviewy(nply
1810           endif           endif
1811        enddo        enddo
1812    
1813   80   continue  c 80   continue
1814          continue
1815                
1816                
1817        if(DEBUG.EQ.1)then        if(DEBUG.EQ.1)then
# Line 2123  cc                                 if(if Line 2139  cc                                 if(if
2139                                         SZX=SZX+ZP(I)*XX                                         SZX=SZX+ZP(I)*XX
2140                                         SSX=SSX+XX                                         SSX=SSX+XX
2141                                         SZ=SZ+ZP(I)                                         SZ=SZ+ZP(I)
2142                                         S1=S1+1.                                                                             S1=S1+1.
2143                                      ENDDO                                      ENDDO
2144                                      DET=SZZ*S1-SZ*SZ                                      DET=SZZ*S1-SZ*SZ
2145                                      AX=(SZX*S1-SZ*SSX)/DET                                      AX=(SZX*S1-SZ*SSX)/DET
# Line 2209  c$$$                                prin Line 2225  c$$$                                prin
2225                       enddo      !end loop on planes  - COPPIA 3                       enddo      !end loop on planes  - COPPIA 3
2226    
2227   31                  continue                       31                  continue                    
2228   1                enddo         !end loop on COPPIA 2  c 1                enddo         !end loop on COPPIA 2
2229                     enddo         !end loop on COPPIA 2
2230                 enddo            !end loop on sensors - COPPIA 2                 enddo            !end loop on sensors - COPPIA 2
2231   20            continue   20            continue
2232              enddo               !end loop on planes  - COPPIA 2              enddo               !end loop on planes  - COPPIA 2
2233    
2234   11         continue  c 11         continue
2235              continue
2236           enddo                  !end loop on COPPIA1           enddo                  !end loop on COPPIA1
2237        enddo                     !end loop on sensors - COPPIA 1        enddo                     !end loop on sensors - COPPIA 1
2238   10   continue   10   continue
# Line 2321  ccccc if(db_used(idbref).eq.1)goto 1188 Line 2339  ccccc if(db_used(idbref).eq.1)goto 1188
2339  *     doublet distance in parameter space  *     doublet distance in parameter space
2340                 distance=                 distance=
2341       $              ((alfayz1(idbref)-alfayz1(idb2))/Dalfayz1)**2       $              ((alfayz1(idbref)-alfayz1(idb2))/Dalfayz1)**2
2342       $              +((alfayz2(idbref)-alfayz2(idb2))/Dalfayz2)**2                     $              +((alfayz2(idbref)-alfayz2(idb2))/Dalfayz2)**2
2343                 distance = sqrt(distance)                 distance = sqrt(distance)
2344                                
2345                 if(distance.lt.cutdistyz)then                 if(distance.lt.cutdistyz)then
# Line 2346  ccccc if(db_used(idbref).eq.1)goto 1188 Line 2364  ccccc if(db_used(idbref).eq.1)goto 1188
2364   1118          continue   1118          continue
2365              enddo               !end loop (2) on DOUBLETS              enddo               !end loop (2) on DOUBLETS
2366                            
2367   1188       continue  c 1188       continue
2368                continue
2369           enddo                  !end loop on... bo?           enddo                  !end loop on... bo?
2370                    
2371           nptloop=npv           nptloop=npv
# Line 2534  c         tr_temp(1)=itr1 Line 2553  c         tr_temp(1)=itr1
2553  *     solo i due parametri spaziali per il momemnto  *     solo i due parametri spaziali per il momemnto
2554                 distance=                 distance=
2555       $              ((alfaxz1(itrref)-alfaxz1(itr2))/Dalfaxz1)**2       $              ((alfaxz1(itrref)-alfaxz1(itr2))/Dalfaxz1)**2
2556       $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2                     $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2
2557                 distance = sqrt(distance)                 distance = sqrt(distance)
2558    
2559    
# Line 2573  c         tr_temp(1)=itr1 Line 2592  c         tr_temp(1)=itr1
2592  11188          continue  11188          continue
2593              enddo               !end loop (2) on TRIPLETS              enddo               !end loop (2) on TRIPLETS
2594                                                
2595  11888       continue  c11888       continue
2596                continue
2597           enddo                  !end loop on... bo?               enddo                  !end loop on... bo?    
2598                    
2599           nptloop=npv           nptloop=npv
# Line 2638  c               mask_view(iv) = 6 Line 2658  c               mask_view(iv) = 6
2658           npt_tot=npt_tot+npt           npt_tot=npt_tot+npt
2659                    
2660           if(DEBUG.EQ.1)then           if(DEBUG.EQ.1)then
2661              print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'                            print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'
2662              print*,'- alfaxz1  ',alfaxz1_av(nclouds_xz)              print*,'- alfaxz1  ',alfaxz1_av(nclouds_xz)
2663              print*,'- alfaxz2  ',alfaxz2_av(nclouds_xz)              print*,'- alfaxz2  ',alfaxz2_av(nclouds_xz)
2664              print*,'- alfaxz3  ',alfaxz3_av(nclouds_xz)              print*,'- alfaxz3  ',alfaxz3_av(nclouds_xz)
# Line 2941  c     $                                 Line 2961  c     $                                
2961                                         xm(nplanes-ip+1)=xPAM                                         xm(nplanes-ip+1)=xPAM
2962                                         ym(nplanes-ip+1)=yPAM                                         ym(nplanes-ip+1)=yPAM
2963                                         zm(nplanes-ip+1)=zPAM                                         zm(nplanes-ip+1)=zPAM
2964                                         resx(nplanes-ip+1)=resxPAM                                                                               resx(nplanes-ip+1)=resxPAM
2965                                         resy(nplanes-ip+1)=resyPAM                                         resy(nplanes-ip+1)=resyPAM
2966                                         if(DEBUG.EQ.1)print*,'(X,Y)'                                         if(DEBUG.EQ.1)print*,'(X,Y)'
2967       $                                      ,nplanes-ip+1,xPAM,yPAM       $                                      ,nplanes-ip+1,xPAM,yPAM
# Line 2952  c     $                                 Line 2972  c     $                                
2972                                         ym_B(nplanes-ip+1) = yPAM_B                                         ym_B(nplanes-ip+1) = yPAM_B
2973                                         zm(nplanes-ip+1)                                         zm(nplanes-ip+1)
2974       $                                      = (zPAM_A+zPAM_B)/2.       $                                      = (zPAM_A+zPAM_B)/2.
2975                                         resx(nplanes-ip+1) = resxPAM                                                                               resx(nplanes-ip+1) = resxPAM
2976                                         resy(nplanes-ip+1) = resyPAM                                         resy(nplanes-ip+1) = resyPAM
2977                                         if(icx.eq.0.and.icy.gt.0)then                                         if(icx.eq.0.and.icy.gt.0)then
2978                                            xgood(nplanes-ip+1)=0.                                            xgood(nplanes-ip+1)=0.
# Line 3011  c                                 chi2=- Line 3031  c                                 chi2=-
3031  *     **********************************************************  *     **********************************************************
3032    
3033                                if(chi2.le.0.)goto 666                                              if(chi2.le.0.)goto 666              
3034                                  if(chi2.ge.1.e08)goto 666 !OPTIMIZATION
3035                                  if(chi2.ne.chi2)goto 666  !OPTIMIZATION
3036    
3037  *     --------------------------  *     --------------------------
3038  *     STORE candidate TRACK INFO  *     STORE candidate TRACK INFO
# Line 3037  c                                    mas Line 3059  c                                    mas
3059    
3060                                   XV_STORE(ip,ntracks)=sngl(xv(ip))                                   XV_STORE(ip,ntracks)=sngl(xv(ip))
3061                                   YV_STORE(ip,ntracks)=sngl(yv(ip))                                   YV_STORE(ip,ntracks)=sngl(yv(ip))
3062                                   ZV_STORE(ip,ntracks)=sngl(zv(ip))                                                                       ZV_STORE(ip,ntracks)=sngl(zv(ip))
3063                                   XM_STORE(ip,ntracks)=sngl(xm(ip))                                   XM_STORE(ip,ntracks)=sngl(xm(ip))
3064                                   YM_STORE(ip,ntracks)=sngl(ym(ip))                                   YM_STORE(ip,ntracks)=sngl(ym(ip))
3065                                   ZM_STORE(ip,ntracks)=sngl(zm(ip))                                   ZM_STORE(ip,ntracks)=sngl(zm(ip))
# Line 3872  c$$$         angy    = 180.*atan(tgtemp) Line 3894  c$$$         angy    = 180.*atan(tgtemp)
3894           if(ip.eq.6.and.ssensor.ne.0)ssensor = 3 - ssensor !notazione paolo x align           if(ip.eq.6.and.ssensor.ne.0)ssensor = 3 - ssensor !notazione paolo x align
3895           LS(IP,ntr)      = ssensor+10*sladder           LS(IP,ntr)      = ssensor+10*sladder
3896    
3897           if(id.ne.0)then  c         if(id.ne.0)then
3898    CCCCCC(10 novembre 2009) PATCH X NUCLEI
3899    C     non ho capito perche', ma durante il ritracciamento dei nuclei
3900    C     (quando una traccia viene trovata ma non e' buona) c'e' qualche variabile
3901    C     che non viene reinizializzata correttamente e i cluster esclusi
3902    C     dal fit risultano ancora inclusi...
3903    C
3904             cltrx(ip,ntr)   = 0
3905             cltry(ip,ntr)   = 0
3906    c$$$         if(
3907    c$$$     $        xgood_nt(ip,ntr).eq.1.and.ygood_nt(ip,ntr).eq.1
3908    c$$$     $        .and.
3909    c$$$     $        id.ne.0)then
3910             if(id.ne.0)then        !patch 30/12/09
3911    
3912  c           >>> is a couple  c           >>> is a couple
3913              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))
3914              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))
# Line 3881  c           >>> is a couple Line 3917  c           >>> is a couple
3917    
3918                 cl_used(cltrx(ip,ntr)) = 1 !tag used clusters                           cl_used(cltrx(ip,ntr)) = 1 !tag used clusters          
3919    
3920                 xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))                               xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))
3921    
3922                 if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0)                 if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0)
3923       $              dedx_x(ip,ntr)=-dedx_x(ip,ntr)       $              dedx_x(ip,ntr)=-dedx_x(ip,ntr)
# Line 3911  c           >>> is a couple Line 3947  c           >>> is a couple
3947                 ypu(ip,ntr)      = corr                 ypu(ip,ntr)      = corr
3948              endif              endif
3949    
3950           elseif(icl.ne.0)then  c$$$         elseif(icl.ne.0)then
3951             endif                  !patch 30/12/09
3952             if(icl.ne.0)then       !patch 30/12/09
3953    
3954              cl_used(icl) = 1    !tag used clusters                        cl_used(icl) = 1    !tag used clusters          
3955    
# Line 4102  c                  call xyz_PAM(0,icl,is Line 4140  c                  call xyz_PAM(0,icl,is
4140                 alfayz2_av_nt(iyz)=alfayz2_av(iyz)                 alfayz2_av_nt(iyz)=alfayz2_av(iyz)
4141                 nnn=nnn+ptcloud_yz(iyz)                 nnn=nnn+ptcloud_yz(iyz)
4142              enddo              enddo
4143              do ipt=1,nnn              do ipt=1,min(ndblt_max_nt,nnn)
4144                 db_cloud_nt(ipt)=db_cloud(ipt)                 db_cloud_nt(ipt)=db_cloud(ipt)
4145               enddo               enddo
4146           endif           endif
# Line 4115  c                  call xyz_PAM(0,icl,is Line 4153  c                  call xyz_PAM(0,icl,is
4153                 alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)                 alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)
4154                 nnn=nnn+ptcloud_xz(ixz)                               nnn=nnn+ptcloud_xz(ixz)              
4155              enddo              enddo
4156              do ipt=1,nnn              do ipt=1,min(ntrpt_max_nt,nnn)
4157                tr_cloud_nt(ipt)=tr_cloud(ipt)                tr_cloud_nt(ipt)=tr_cloud(ipt)
4158               enddo               enddo
4159           endif           endif
Legend:
Removed from v.1.37  
changed lines
  Added in v.1.44
  ViewVC Help
Powered by ViewVC 1.1.23