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

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

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

revision 1.1 by mocchiut, Fri May 19 13:15:55 2006 UTC revision 1.10 by mocchiut, Thu Jan 16 15:29:52 2014 UTC
# Line 7  Line 7 
7  *            *          
8  *************************************************************************  *************************************************************************
9    
10        subroutine cncomp(i,j)    !(view, VA1)        subroutine cncomp(i,j,errflag)    !(view, VA1)
11    
12        include 'commontracker.f'        include 'commontracker.f'
13          include 'level1.f'
14        include 'common_reduction.f'        include 'common_reduction.f'
15        include 'calib.f'        include 'calib.f'
16          
17    
18        integer errflag           !error flag to mark no signal free VA1        integer errflag           !error flag to mark no signal free VA1
   
19        integer clstr_old(nstrips_va1) !flag storage vector        integer clstr_old(nstrips_va1) !flag storage vector
20    c     EM GCC4.7 - trying to avoind implicit conversion from REAL(8) to REAL(4) so use double everywhere in routines and cast to float only when
21        real signal(nstrips_va1)  !"signal" (=adc-ped) value storage vector  c     variables from common includes are used
22          real(8) signal(nstrips_va1)  !"signal" (=adc-ped) value storage vector
23        real smean, ssigma        !"signal" mean and sigma        real(8) smean, ssigma        !"signal" mean and sigma
24        real cut                  !"strange" strip exclusion cut        real(8) cut                  !"strange" strip exclusion cut
25    c     END EM
26        integer newclstr          !flag to warn about new found clusters to be        integer newclstr          !flag to warn about new found clusters to be
27                                  ! excluded from common noise computation  c                               ! excluded from common noise computation
   
   
 c     call HBOOK1(20000+100*i+j,' ',30,0.,30.,0.) !???  
28    
29  c------------------------------------------------------------------------  c------------------------------------------------------------------------
30  c      c
31  c     variables initialization  c     variables initialization
32  c      c
33  c------------------------------------------------------------------------  c------------------------------------------------------------------------
34        do k=1,nstrips_va1        !loops on strips        do k=1,nstrips_va1        !loops on strips
35           clstr(i,j,k)=1         !initializes signal affected strips flag           clstr(i,j,k)=1         !initializes signal affected strips flag
# Line 39  c--------------------------------------- Line 37  c---------------------------------------
37           strange(i,j,k)=1       !initializes unusually high or low signal           strange(i,j,k)=1       !initializes unusually high or low signal
38        enddo                     ! affected strips flag        enddo                     ! affected strips flag
39    
40        newclstr=1                !flag to warn about new found signal  c------------------------------------------------------------------------
41                                  ! affected strips  c     (september 2007)
42    c     remove from CN computation the first and the last 3 channels of
43    c     each X view, becouse they ar not connected to any strip
44    c------------------------------------------------------------------------
45          if(mod(i,2).eq.0)then
46             if(j.eq.1)then
47                do k=1,3
48                   strange(i,j,k)=0
49                enddo
50             elseif(j.eq.nva1_ladder)then
51                do k=nstrips_va1,nstrips_va1-2,-1
52                   strange(i,j,k)=0
53                enddo
54             endif
55          endif
56    
57          newclstr=1                !flag to warn about new found signal
58    c                               ! affected strips
59  c------------------------------------------------------------------------  c------------------------------------------------------------------------
60  c      c
61  c     high or low signal affected strips exclusion: computes "signal" (=adc-ped)  c     high or low signal affected strips exclusion: computes "signal" (=adc-ped)
62  c     mean value and sigma, and cuts from common noise computation strips  c     mean value and sigma, and cuts from common noise computation strips
63  c     whose ABS(signal) exceeds scut*sigma  c     whose ABS(signal) exceeds scut*sigma
64  c      c
65  c------------------------------------------------------------------------  c------------------------------------------------------------------------
66        countme=0                 !???        countme=0                 !???
67   666  continue                  !???   666  continue                  !???
# Line 59  c--------------------------------------- Line 71  c---------------------------------------
71        nstr=0        nstr=0
72                
73        do k=1,nstrips_va1        do k=1,nstrips_va1
74           nstr=nstr+strange(i,j,k) !uses only           nstr = nstr + strange(i,j,k) !uses only
75           if(mod(i,2).eq.1) then !odd strip ---> Y view           if(mod(i,2).eq.1) then ! ---> Y view
76              signal(k)= - (DBLE(adc(i,j,k))-pedestal(i,j,k)) !negative signal              signal(k) = - (DBLE(adc(i,j,k)) - pedestal(i,j,k)) !negative signal
77           else                   !even strip ---> X view           else                   ! ---> X view
78              signal(k)= DBLE(adc(i,j,k))-pedestal(i,j,k) !positive signal              signal(k) =    DBLE(adc(i,j,k)) - pedestal(i,j,k) !positive signal
79           endif           endif
80                     smean = smean + signal(k)*strange(i,j,k)
81           smean=smean+signal(k)*strange(i,j,k)           ssigma = ssigma + (signal(k)**2)*strange(i,j,k)
          ssigma=ssigma+(signal(k)**2)*strange(i,j,k)  
           
 c     call HFILL(10000+100*i+j,signal(k),0.,1.) !???  
82        enddo        enddo
83                
84        smean=smean/nstr          !strips value distribution mean        smean=smean/nstr          !strips value distribution mean      
         
85        ssigma=SQRT((ssigma/nstr)-smean**2) !strips value distribution sigma        ssigma=SQRT((ssigma/nstr)-smean**2) !strips value distribution sigma
86                
87        cut=scut*ssigma           !exclusion cut        cut=scut*ssigma           !exclusion cut
88                
89          nco=0
90          nbo=0
91        do k=1,nstrips_va1        do k=1,nstrips_va1
 c     call HFILL(20000+100*i+j,ABS(signal(k)-smean)/ssigma,0.,1.) !???  
92           if(ABS(signal(k)-smean).gt.cut) then           if(ABS(signal(k)-smean).gt.cut) then
 c     print*,i,j,k,signal(k),abs(signal(k)),cut,strange(i,j,k) !???  
93              strange(i,j,k)=0    !marks strips exceeding cut              strange(i,j,k)=0    !marks strips exceeding cut
94    c            print*,i,j,k,signal(k),smean
95           endif           endif
96             nco=nco+strange(i,j,k)
97             nbo=nbo+bad(i,j,k)
98        enddo                     ! in order not to use them in CN computation        enddo                     ! in order not to use them in CN computation
99    
100    c$$$      if(i.eq.12.and.(j.eq.2.or.j.eq.3))then
101    c$$$         print*,'view ',i,' vk ',j
102    c$$$         print*,'ADC (1-51-128) = ',adc(i,j,1),adc(i,j,52),adc(i,j,128)
103    c$$$         print*,'<ADC-PED> = ',smean
104    c$$$         print*,'s         = ',ssigma
105    c$$$         print*,'nstrange  = ',128-nco
106    c$$$         print*,'nbad      = ',128-nbo
107    c$$$      endif
108    
109        countme=countme+1         !???        countme = countme + 1         !???
110        if (countme.le.3) goto 666 !???        if (countme.le.3) goto 666 !???
111    
   
112  c------------------------------------------------------------------------  c------------------------------------------------------------------------
113  c      c    
114  c     common noise computation  c     common noise computation
115  c      c    
116  c------------------------------------------------------------------------  c-----------------------------------------------------------------------
117        do while(newclstr.eq.1)   !loops on this VA1 till no new signal  *     loops on this VA1 till no new signal affected strips are found
118                                  ! affected strips are found        do while(newclstr.eq.1)  
119    
120           newclstr=0             !to exit from loop if no new cluster is           newclstr=0             !to exit from loop if no new cluster is found
                                 ! found  
121                    
122           errflag=0           errflag=0
123             call cnoise(i,j,errflag) !(view, VA1, error flag) computes cn
124           call cnoise(i,j,errflag) !(view, VA1, error flag) computes common           if(errflag.eq.1) goto 10 !goes to next VA1: this one has no signal-free strips...
                                 ! noise  
   
 c     print*,cn(i,j)         !???  
           
          if(errflag.eq.1) goto 10 !goes to next VA1: this one has no signal  
                                 ! free strips...  
125                    
126           call cutcn(i,j)        !(view, VA1) excludes clusters from           call cutcn(i,j)        !(view, VA1) excludes clusters from cn computation
                                 ! common noise calculation  
127                    
128           ncs=0                  !initializes number of strips not excluded by cncut           ncs=0                  !initializes number of strips not excluded by cncut
   
129           do k=1,nstrips_va1     !loops on strips           do k=1,nstrips_va1     !loops on strips
130              if(clstr(i,j,k).ne.clstr_old(k)) then !checks if there are  *           checks if there are new found clusters, and if so sets
131                                  ! new found clusters, and if so sets              if(clstr(i,j,k).ne.clstr_old(k)) then
132                 newclstr=1       ! newclstr flag = 1                 newclstr=1                      
133                                 clstr_old(k)=clstr(i,j,k)  
134                 clstr_old(k)=clstr(i,j,k) !stores cluster flags in              endif              
             endif               ! clstr_old variable  
   
135              iok=strange(i,j,k)*bad(i,j,k)*clstr(i,j,k)              iok=strange(i,j,k)*bad(i,j,k)*clstr(i,j,k)
   
136              ncs=ncs+iok         !counts number of good strips for cn computation              ncs=ncs+iok         !counts number of good strips for cn computation
   
137           enddo           enddo
138    
139        enddo                     !ends do while loop when there are no new        enddo                     !ends do while
                                 ! clusters  
   
 c      call HFILL(666,FLOAT(ncs),0.,1.) !???  
   
   
 c$$$      if(ncs.lt.20) then        !warns if too many strips have been excluded from CN  
 c$$$                                ! computation  
 c$$$         print*,'cncomp: WARNING, LESS THAN 20 STRIPS PASSED CN CUT'  
 c$$$     $        //' ON VA1 ',j,', VIEW ',i !NB questo errore e' "un po'" in conflitto  
 c$$$                                ! con quello che setta errflag (vedi cnoise.f)...  
 c$$$  
 c$$$      endif  
140    
141   10   continue   10   continue
142    
# Line 174  c$$$      endif Line 167  c$$$      endif
167        subroutine cnoise(i,j,gulp) !(view, VA1)        subroutine cnoise(i,j,gulp) !(view, VA1)
168    
169        include 'commontracker.f'        include 'commontracker.f'
170          include 'level0.f'
171          include 'level1.f'
172        include 'common_reduction.f'        include 'common_reduction.f'
173        include 'calib.f'        include 'calib.f'
174                
# Line 182  c$$$      endif Line 177  c$$$      endif
177                
178        ncn=0                     !number of strips in cn computation        ncn=0                     !number of strips in cn computation
179        cn(i,j)=0                 !initializes cn variable        cn(i,j)=0                 !initializes cn variable
180          cnrms(i,j)=0              !initializes cn rms
181          cnn(i,j)=0                !initialize cn flag
182    
183        do k=1,nstrips_va1        !loops on strips        do k=1,nstrips_va1        !loops on strips
184           iok=strange(i,j,k)*bad(i,j,k)*clstr(i,j,k) !flag to mark strange, bad  *        tags strange, bad or signal-affected strips
185                                  ! or signal affected strips           iok = strange(i,j,k)*bad(i,j,k)*clstr(i,j,k)
186  ccc         print*,i,j,k,strange(i,j,k),bad(i,j,k),clstr(i,j,k),iok !???           cn(i,j) = cn(i,j)  
187         $        + REAL((adc(i,j,k)) - pedestal(i,j,k))*REAL(iok) ! EM GCC4.7 CN IS DEFINED AS REAL NOT DOUBLE
188           cn(i,j)=cn(i,j) + (DBLE(adc(i,j,k)) - pedestal(i,j,k))*iok !sums ADC-PED           cnrms(i,j) = cnrms(i,j)  !EM GCC4.7 CNRMS IS DEFINED AS REAL NOT DOUBLE
189                                  ! values to compute common noise       $        + (REAL(adc(i,j,k)) - pedestal(i,j,k)) !EM GCC4.7 CNRMS IS DEFINED AS REAL NOT DOUBLE
190         $        *(REAL(adc(i,j,k)) - pedestal(i,j,k))*REAL(iok) !EM GCC4.7 CNRMS IS DEFINED AS REAL NOT DOUBLE
191           ncn = ncn + iok            !counts number of strips in cn computation           ncn = ncn + iok            !counts number of strips in cn computation
192        enddo        enddo
193          
194  ccc      print*,'ncn= ',ncn        if(ncn.lt.NSTRIPMIN) then         !no signal free strips on this VA1...
195        if(ncn.eq.0) then         !no signal free strips on this VA1...           if(ncn.eq.0)then
196           print*,'cnoise: WARNING, NO SIGNAL FREE STRIPS ON VA1 ',j,              if(debug.eq.1)print*,' WARNING - cnoise: ',
197       $        ', VIEW ',i       $        'no strips for CN computation on VA1 ',j,
198         $        ', VIEW ',i,'  >>> FAILED '
199             else
200                if(debug.eq.1)print*,' WARNING - cnoise: ',
201         $        'less than ',NSTRIPMIN
202         $           ,' strips for CN computation on VA1 ',j,
203         $        ', VIEW ',i,'  >>> FAILED '
204             endif
205           gulp=1           gulp=1
206             cnn(i,j) = 0
207        else        else
208           cn(i,j)=cn(i,j)/DBLE(ncn) !computes common noise           cn(i,j)=cn(i,j)/REAL(ncn) !<<<< computes common noise  EM GCC4.7 CN IS REAL NOT DOUBLE
209           gulp=0                 !resets error flag           cnrms(i,j)= SQRT( cnrms(i,j)/REAL(ncn) - cn(i,j)**2 ) ! EM GCC4.7 CN IS REAL NOT DOUBLE
210             cnn(i,j) = ncn
211             gulp=0                
212    c$$$         print*,'Event ',eventn(1)
213    c$$$     $        ,': cn(',i,',',j,')= ',cn(i,j),' ncn ',ncn
214            
215             if(debug.eq.1.and.ABS(cn(i,j)).gt.1000)
216         $        print*,'Event ',eventn(1)
217         $        ,': cn(',i,',',j,')= ',cn(i,j),' ncn ',ncn
218        endif        endif
219    
220        return        return
# Line 227  ccc      print*,'ncn= ',ncn Line 241  ccc      print*,'ncn= ',ncn
241        subroutine cutcn(i,j)     !(view, VA1)        subroutine cutcn(i,j)     !(view, VA1)
242    
243        include 'commontracker.f'        include 'commontracker.f'
244          include 'level1.f'
245        include 'common_reduction.f'        include 'common_reduction.f'
246        include 'calib.f'        include 'calib.f'
247    
# Line 238  ccc      print*,'ncn= ',ncn Line 253  ccc      print*,'ncn= ',ncn
253                                  ! seed                                  ! seed
254        integer ir, il            !flags to exit loop on reaching VA1 extremes        integer ir, il            !flags to exit loop on reaching VA1 extremes
255    
256        real valuec                !cluster seed signal        real(8) valuec                !cluster seed signal   EM GCC4.7
257        real cut,stripcut         !cluster seed cut        real cut,stripcut         !cluster seed cut
258    
259        real valuel, valuer       !left and right strips signal        real(8) valuel, valuer       !left and right strips signal   EM GCC4.7
260        real stripcnincut         !strip include cut        real stripcnincut         !strip include cut
261    
262        skip = 0                  !initializes skip        skip = 0                  !initializes skip
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.23