/[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.41 by mocchiut, Tue Aug 11 12:58:24 2009 UTC revision 1.46 by mocchiut, Fri Jan 17 10:14:38 2014 UTC
# Line 528  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 545  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 1369  c--------------------------------------- Line 1379  c---------------------------------------
1379                 iv1=iside                 iv1=iside
1380                 iv2=mod(iside,4)+1                 iv2=mod(iside,4)+1
1381  *     straight line passing trhough two consecutive vertexes  *     straight line passing trhough two consecutive vertexes
1382                 AA = (yvv(iv1)-yvv(iv2))/(xvv(iv1)-xvv(iv2))                 AA = REAL((yvv(iv1)-yvv(iv2))/(xvv(iv1)-xvv(iv2))) !EM GCC4.7
1383                 BB = yvv(iv1) - AA*xvv(iv1)                 BB = REAL(yvv(iv1) - AA*xvv(iv1)) !EM GCC4.7
1384  *     point along the straight line closer to the track  *     point along the straight line closer to the track
1385                 xoo = (xPAM+AA*yPAM-AA*BB)/(1+AA**2)                 xoo = (xPAM+AA*yPAM-AA*BB)/(1+AA**2)
1386                 yoo = AA*xoo + BB                 yoo = AA*xoo + BB
# Line 1382  c--------------------------------------- Line 1392  c---------------------------------------
1392                 iv1=iside                 iv1=iside
1393                 iv2=mod(iside,4)+1                 iv2=mod(iside,4)+1
1394  *     straight line passing trhough two consecutive vertexes  *     straight line passing trhough two consecutive vertexes
1395                 AA = (xvv(iv1)-xvv(iv2))/(yvv(iv1)-yvv(iv2))                 AA = REAL((xvv(iv1)-xvv(iv2))/(yvv(iv1)-yvv(iv2))) !EM GCC4.7
1396                 BB = xvv(iv1) - AA*yvv(iv1)                 BB = REAL(xvv(iv1) - AA*yvv(iv1)) !EM GCC4.7
1397  *     point along the straight line closer to the track  *     point along the straight line closer to the track
1398                 yoo = (yPAM+AA*xPAM-AA*BB)/(1+AA**2)                 yoo = (yPAM+AA*xPAM-AA*BB)/(1+AA**2)
1399                 xoo = AA*yoo + BB                 xoo = AA*yoo + BB
# Line 1966  c$$$               print*,'(1) ip ',ip1, Line 1976  c$$$               print*,'(1) ip ',ip1,
1976  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)
1977  c               call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)  c               call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)
1978                 call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.)                 call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.)
1979                 xm1=xPAM                 xm1=REAL(xPAM) !EM GCC4.7
1980                 ym1=yPAM                 ym1=REAL(yPAM) !EM GCC4.7
1981                 zm1=zPAM                                   zm1=REAL(zPAM) !EM GCC4.7
1982    
1983                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
1984  c$$$                  print*,'(2) ip ',ip2  c$$$                  print*,'(2) ip ',ip2
# Line 1989  c                        call xyz_PAM Line 1999  c                        call xyz_PAM
1999  c     $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)  c     $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)
2000                          call xyz_PAM                          call xyz_PAM
2001       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.)       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.)
2002                          xm2=xPAM                          xm2=REAL(xPAM) !EM GCC4.7
2003                          ym2=yPAM                          ym2=REAL(yPAM) !EM GCC4.7
2004                          zm2=zPAM                          zm2=REAL(zPAM) !EM GCC4.7
2005    
2006  *                       ---------------------------------------------------  *                       ---------------------------------------------------
2007  *                       both couples must have a y-cluster  *                       both couples must have a y-cluster
# Line 2032  ccc                        print*,'<doub Line 2042  ccc                        print*,'<doub
2042  *     tg(th_yz)  *     tg(th_yz)
2043                          alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2)                          alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2)
2044  *     y0 (cm)  *     y0 (cm)
2045                          alfayz1(ndblt)=alfayz2(ndblt)*(zini-zm1)+ym1                      alfayz1(ndblt)=alfayz2(ndblt)*(REAL(zini)-zm1)+ym1! EM GCC4.7 zm1, ym1 and alfayz1/2 are REAL
2046                                                    
2047  ****  -----------------------------------------------****  ****  -----------------------------------------------****
2048  ****  reject non phisical couples                    ****  ****  reject non phisical couples                    ****
# Line 2095  c     $                               (i Line 2105  c     $                               (i
2105                                   call xyz_PAM                                   call xyz_PAM
2106       $                                (icx3,icy3,is3,PFAdef,PFAdef       $                                (icx3,icy3,is3,PFAdef,PFAdef
2107       $                                ,0.,0.,0.,0.)       $                                ,0.,0.,0.,0.)
2108                                   xm3=xPAM                                   xm3=REAL(xPAM) !EM GCC4.7
2109                                   ym3=yPAM                                   ym3=REAL(yPAM) !EM GCC4.7
2110                                   zm3=zPAM                                   zm3=REAL(zPAM) !EM GCC4.7
2111    
2112    
2113  *     find the circle passing through the three points  *     find the circle passing through the three points
# Line 2134  cc                                 if(if Line 2144  cc                                 if(if
2144                                      DET=SZZ*S1-SZ*SZ                                      DET=SZZ*S1-SZ*SZ
2145                                      AX=(SZX*S1-SZ*SSX)/DET                                      AX=(SZX*S1-SZ*SSX)/DET
2146                                      BX=(SZZ*SSX-SZX*SZ)/DET                                      BX=(SZZ*SSX-SZX*SZ)/DET
2147                                      X0  = AX*ZINI+BX                                      X0  = AX*REAL(ZINI)+BX ! EM GCC4.7
2148                                                                            
2149                                   endif                                   endif
2150    
# Line 2173  ccc                                 prin Line 2183  ccc                                 prin
2183                                                                    
2184                                   if(radius.ne.0.and.xc.lt.0)then                                   if(radius.ne.0.and.xc.lt.0)then
2185  *************POSITIVE DEFLECTION  *************POSITIVE DEFLECTION
2186               alfaxz1(ntrpt) = xc+sqrt(radius**2-(ZINI-zc)**2)             alfaxz1(ntrpt) = xc+sqrt(radius**2-(REAL(ZINI)-zc)**2) !EM GCC4.7
2187               alfaxz2(ntrpt) = (ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)             alfaxz2(ntrpt) = (REAL(ZINI)-zc)/
2188               alfaxz3(ntrpt) = 1/radius       $          sqrt(radius**2-(REAL(ZINI)-zc)**2) !EM GCC4.7
2189               alfaxz3(ntrpt) = 1/radius
2190                                  else if(radius.ne.0.and.xc.ge.0)then                                  else if(radius.ne.0.and.xc.ge.0)then
2191  *************NEGATIVE DEFLECTION  *************NEGATIVE DEFLECTION
2192               alfaxz1(ntrpt) = xc-sqrt(radius**2-(ZINI-zc)**2)             alfaxz1(ntrpt) = xc-sqrt(radius**2-(REAL(ZINI)-zc)**2)
2193               alfaxz2(ntrpt) = -(ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)             alfaxz2(ntrpt) = -(REAL(ZINI)-zc)/
2194               alfaxz3(ntrpt) = -1/radius       $          sqrt(radius**2-(REAL(ZINI)-zc)**2) !EM GCC4.7
2195               alfaxz3(ntrpt) = -1/radius
2196                                  else if(radius.eq.0)then                                  else if(radius.eq.0)then
2197  *************straight fit  *************straight fit
2198               alfaxz1(ntrpt) = X0               alfaxz1(ntrpt) = X0
# Line 3094  c                                    mas Line 3106  c                                    mas
3106                                   enddo                                   enddo
3107                                enddo                                enddo
3108                                                                
3109                                RCHI2_STORE(ntracks)=chi2                                RCHI2_STORE(ntracks)=REAL(chi2)
3110                                                                
3111  *     --------------------------------  *     --------------------------------
3112  *     STORE candidate TRACK INFO - end  *     STORE candidate TRACK INFO - end
# Line 3162  cc         return Line 3174  cc         return
3174        character*10 PFA        character*10 PFA
3175        common/FINALPFA/PFA        common/FINALPFA/PFA
3176    
3177          double precision xmm,rxmm,xmm_A,xmm_B !EM GCC4.7
3178          double precision ymm,rymm,ymm_A,ymm_B !EM GCC4.7
3179          double precision zmm,zmm_A,zmm_B !EM GCC4.7
3180          double precision clincnewc !EM GCC4.7
3181          double precision clincnew !EM GCC4.7
3182    
3183        real k(6)        real k(6)
3184        DATA k/1.099730,0.418900,0.220939,0.220907,0.418771,1.100674/        DATA k/1.099730,0.418900,0.220939,0.220907,0.418771,1.100674/
3185    
# Line 3380  c               distance = distance / RC Line 3398  c               distance = distance / RC
3398                    idm = id                                      idm = id                  
3399                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3400                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3401                    clincnewc=10*sqrt(rymm**2+rxmm**2                    clincnewc=10.*dsqrt(rymm**2+rxmm**2
3402       $                 +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2)))       $            +DBLE(RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2))))! EM GCC4.7
3403                 endif                 endif
3404   1188          continue   1188          continue
3405              enddo               !end loop on couples on plane icp              enddo               !end loop on couples on plane icp
# Line 3549  c               distance = distance / RC Line 3567  c               distance = distance / RC
3567              if(iclm.ne.0)then              if(iclm.ne.0)then
3568                 if(mod(VIEW(iclm),2).eq.0)then                 if(mod(VIEW(iclm),2).eq.0)then
3569                    clincnew=                    clincnew=
3570       $                 20*       $            20.*     !EM GCC4.7
3571       $                 sqrt(rxmm**2+RCHI2_STORE(ibest)*k(ip)*cov(1,1))       $            dsqrt(rxmm**2 +
3572         $             DBLE(RCHI2_STORE(ibest)*k(ip))*cov(1,1))
3573                 else if(mod(VIEW(iclm),2).ne.0)then                 else if(mod(VIEW(iclm),2).ne.0)then
3574                    clincnew=                    clincnew=
3575       $                 10*       $            10.* !EM GCC4.7
3576       $                 sqrt(rymm**2+RCHI2_STORE(ibest)*k(ip)*cov(2,2))       $            dsqrt(rymm**2 +
3577         $             DBLE(RCHI2_STORE(ibest)*k(ip))*cov(2,2))
3578                 endif                 endif
3579    
3580                 if(distmin.le.clincnew)then                   if(distmin.le.clincnew)then  
# Line 3803  c               distance = distance / RC Line 3823  c               distance = distance / RC
3823        character*10 PFA        character*10 PFA
3824        common/FINALPFA/PFA        common/FINALPFA/PFA
3825    
3826        real sinth,phi,pig        real sinth,phi,pig, npig ! EM GCC4.7
3827        integer ssensor,sladder        integer ssensor,sladder
3828        pig=acos(-1.)        pig=acos(-1.)
3829    
# Line 3813  c               distance = distance / RC Line 3833  c               distance = distance / RC
3833        chi2_nt(ntr)        = sngl(chi2)        chi2_nt(ntr)        = sngl(chi2)
3834        nstep_nt(ntr)       = nstep        nstep_nt(ntr)       = nstep
3835  *     -------------------------------------  *     -------------------------------------
3836        phi   = al(4)                  phi   = REAL(al(4))
3837        sinth = al(3)                    sinth = REAL(al(3))
3838        if(sinth.lt.0)then              if(sinth.lt.0)then      
3839           sinth = -sinth                   sinth = -sinth        
3840           phi = phi + pig                 phi = phi + pig      
# Line 3884  c$$$         angy    = 180.*atan(tgtemp) Line 3904  c$$$         angy    = 180.*atan(tgtemp)
3904           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
3905           LS(IP,ntr)      = ssensor+10*sladder           LS(IP,ntr)      = ssensor+10*sladder
3906    
3907           if(id.ne.0)then  c         if(id.ne.0)then
3908    CCCCCC(10 novembre 2009) PATCH X NUCLEI
3909    C     non ho capito perche', ma durante il ritracciamento dei nuclei
3910    C     (quando una traccia viene trovata ma non e' buona) c'e' qualche variabile
3911    C     che non viene reinizializzata correttamente e i cluster esclusi
3912    C     dal fit risultano ancora inclusi...
3913    C
3914             cltrx(ip,ntr)   = 0
3915             cltry(ip,ntr)   = 0
3916    c$$$         if(
3917    c$$$     $        xgood_nt(ip,ntr).eq.1.and.ygood_nt(ip,ntr).eq.1
3918    c$$$     $        .and.
3919    c$$$     $        id.ne.0)then
3920             if(id.ne.0)then        !patch 30/12/09
3921    
3922  c           >>> is a couple  c           >>> is a couple
3923              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))
3924              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))
# Line 3923  c           >>> is a couple Line 3957  c           >>> is a couple
3957                 ypu(ip,ntr)      = corr                 ypu(ip,ntr)      = corr
3958              endif              endif
3959    
3960           elseif(icl.ne.0)then  c$$$         elseif(icl.ne.0)then
3961             endif                  !patch 30/12/09
3962             if(icl.ne.0)then       !patch 30/12/09
3963    
3964              cl_used(icl) = 1    !tag used clusters                        cl_used(icl) = 1    !tag used clusters          
3965    
# Line 4021  c         if( mask_view(iv).ne.0 )good2( Line 4057  c         if( mask_view(iv).ne.0 )good2(
4057                 do is=1,2                 do is=1,2
4058  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
4059  c                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)  c                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
4060                    call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.,0.,0.)                    call xyz_PAM(icl,0,is,PFAdef,'    ',0.,0.,0.,0.)
4061                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2                    xs(is,nclsx) = REAL((xPAM_A+xPAM_B)/2.) ! EM GCC4.7
4062                 enddo                 enddo
4063              else                          !=== Y views              else                          !=== Y views
4064                 nclsy = nclsy + 1                 nclsy = nclsy + 1
# Line 4037  c                  call xyz_PAM(icl,0,is Line 4073  c                  call xyz_PAM(icl,0,is
4073                 do is=1,2                 do is=1,2
4074  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
4075  c                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)  c                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
4076                    call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.,0.,0.)                    call xyz_PAM(0,icl,is,'    ',PFAdef,0.,0.,0.,0.)
4077                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2                    ys(is,nclsy) = REAL((yPAM_A+yPAM_B)/2.) ! EM GCC4.7
4078                 enddo                 enddo
4079              endif              endif
4080           endif           endif
Legend:
Removed from v.1.41  
changed lines
  Added in v.1.46
  ViewVC Help
Powered by ViewVC 1.1.23