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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.10 - (hide annotations) (download)
Thu Jan 16 15:29:52 2014 UTC (10 years, 10 months ago) by mocchiut
Branch: MAIN
CVS Tags: v10RED, v10REDr01, HEAD
Changes since 1.9: +16 -11 lines
Compilation warnings using GCC4.7 fixed

1 mocchiut 1.1 *************************************************************************
2     *
3     * Subroutine cncomp.f
4     *
5     * iterates common noise computation subroutine (./cnoise.f) and cluster
6     * cutting subroutine (./cutcn.f) till no more clusters are found
7     *
8     *************************************************************************
9    
10 pam-fi 1.2 subroutine cncomp(i,j,errflag) !(view, VA1)
11 mocchiut 1.1
12     include 'commontracker.f'
13 pam-fi 1.4 include 'level1.f'
14 mocchiut 1.1 include 'common_reduction.f'
15     include 'calib.f'
16 mocchiut 1.10
17 mocchiut 1.1
18     integer errflag !error flag to mark no signal free VA1
19     integer clstr_old(nstrips_va1) !flag storage vector
20 mocchiut 1.10 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     c variables from common includes are used
22     real(8) signal(nstrips_va1) !"signal" (=adc-ped) value storage vector
23     real(8) smean, ssigma !"signal" mean and sigma
24     real(8) cut !"strange" strip exclusion cut
25     c END EM
26 mocchiut 1.1 integer newclstr !flag to warn about new found clusters to be
27 pam-fi 1.2 c ! excluded from common noise computation
28 mocchiut 1.1
29     c------------------------------------------------------------------------
30 pam-fi 1.2 c
31 mocchiut 1.1 c variables initialization
32 pam-fi 1.2 c
33 mocchiut 1.1 c------------------------------------------------------------------------
34     do k=1,nstrips_va1 !loops on strips
35     clstr(i,j,k)=1 !initializes signal affected strips flag
36     clstr_old(k)=1 !initializes signal affected strips storage
37     strange(i,j,k)=1 !initializes unusually high or low signal
38     enddo ! affected strips flag
39    
40 pam-fi 1.9 c------------------------------------------------------------------------
41     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 mocchiut 1.1 newclstr=1 !flag to warn about new found signal
58 pam-fi 1.2 c ! affected strips
59 mocchiut 1.1 c------------------------------------------------------------------------
60 pam-fi 1.2 c
61 mocchiut 1.1 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
63     c whose ABS(signal) exceeds scut*sigma
64 pam-fi 1.2 c
65 mocchiut 1.1 c------------------------------------------------------------------------
66     countme=0 !???
67     666 continue !???
68    
69     smean=0. !initialization
70     ssigma=0.
71     nstr=0
72    
73     do k=1,nstrips_va1
74 pam-fi 1.2 nstr = nstr + strange(i,j,k) !uses only
75 pam-fi 1.6 if(mod(i,2).eq.1) then ! ---> Y view
76 pam-fi 1.2 signal(k) = - (DBLE(adc(i,j,k)) - pedestal(i,j,k)) !negative signal
77 pam-fi 1.6 else ! ---> X view
78 pam-fi 1.2 signal(k) = DBLE(adc(i,j,k)) - pedestal(i,j,k) !positive signal
79 mocchiut 1.1 endif
80 pam-fi 1.2 smean = smean + signal(k)*strange(i,j,k)
81     ssigma = ssigma + (signal(k)**2)*strange(i,j,k)
82 mocchiut 1.1 enddo
83    
84 pam-fi 1.6 smean=smean/nstr !strips value distribution mean
85 mocchiut 1.1 ssigma=SQRT((ssigma/nstr)-smean**2) !strips value distribution sigma
86    
87     cut=scut*ssigma !exclusion cut
88    
89 pam-fi 1.3 nco=0
90     nbo=0
91 mocchiut 1.1 do k=1,nstrips_va1
92     if(ABS(signal(k)-smean).gt.cut) then
93     strange(i,j,k)=0 !marks strips exceeding cut
94 pam-fi 1.3 c print*,i,j,k,signal(k),smean
95 mocchiut 1.1 endif
96 pam-fi 1.3 nco=nco+strange(i,j,k)
97     nbo=nbo+bad(i,j,k)
98 mocchiut 1.1 enddo ! in order not to use them in CN computation
99    
100 pam-fi 1.3 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 mocchiut 1.1
109 pam-fi 1.2 countme = countme + 1 !???
110 mocchiut 1.1 if (countme.le.3) goto 666 !???
111    
112     c------------------------------------------------------------------------
113     c
114     c common noise computation
115     c
116 pam-fi 1.2 c-----------------------------------------------------------------------
117     * loops on this VA1 till no new signal affected strips are found
118     do while(newclstr.eq.1)
119 mocchiut 1.1
120 pam-fi 1.2 newclstr=0 !to exit from loop if no new cluster is found
121 mocchiut 1.1
122     errflag=0
123 pam-fi 1.2 call cnoise(i,j,errflag) !(view, VA1, error flag) computes cn
124     if(errflag.eq.1) goto 10 !goes to next VA1: this one has no signal-free strips...
125 mocchiut 1.1
126 pam-fi 1.2 call cutcn(i,j) !(view, VA1) excludes clusters from cn computation
127 mocchiut 1.1
128     ncs=0 !initializes number of strips not excluded by cncut
129     do k=1,nstrips_va1 !loops on strips
130 pam-fi 1.2 * checks if there are new found clusters, and if so sets
131     if(clstr(i,j,k).ne.clstr_old(k)) then
132     newclstr=1
133     clstr_old(k)=clstr(i,j,k)
134     endif
135 mocchiut 1.1 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
137     enddo
138    
139 pam-fi 1.2 enddo !ends do while
140 mocchiut 1.1
141     10 continue
142    
143     return
144     end
145    
146    
147    
148    
149     *************************************************************************
150     *
151     * Subroutine cnoise.f!DA COMMENTARE!???
152     *
153     * uses adc(nviews,nva1_view,nstrips_va1) and
154     * pedestal(nviews,nva1_view,nstrips_va1) variables to compute common noise,
155     * and fills cn(nviews,nva1_view) variable. in the computation only
156     * not-bad and not-signal-affected strips are used
157     * (bad(nviews,nva1_view,nstrips_va1) and
158     * clstr(nviews,nva1_view,nstrips_va1) flags)
159     *
160     * needs:
161     * - ./common_calib.f
162     *
163     * to be called inside ./cncomp.f
164     *
165     *************************************************************************
166    
167     subroutine cnoise(i,j,gulp) !(view, VA1)
168    
169     include 'commontracker.f'
170 pam-fi 1.4 include 'level0.f'
171     include 'level1.f'
172 mocchiut 1.1 include 'common_reduction.f'
173     include 'calib.f'
174    
175    
176     integer gulp !error flag
177    
178     ncn=0 !number of strips in cn computation
179     cn(i,j)=0 !initializes cn variable
180 pam-fi 1.6 cnrms(i,j)=0 !initializes cn rms
181 pam-fi 1.4 cnn(i,j)=0 !initialize cn flag
182 mocchiut 1.1
183     do k=1,nstrips_va1 !loops on strips
184 pam-fi 1.2 * tags strange, bad or signal-affected strips
185 pam-fi 1.3 iok = strange(i,j,k)*bad(i,j,k)*clstr(i,j,k)
186 mocchiut 1.10 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     cnrms(i,j) = cnrms(i,j) !EM GCC4.7 CNRMS IS DEFINED AS REAL NOT DOUBLE
189     $ + (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 mocchiut 1.1 ncn = ncn + iok !counts number of strips in cn computation
192     enddo
193 pam-fi 1.2
194 pam-fi 1.3 if(ncn.lt.NSTRIPMIN) then !no signal free strips on this VA1...
195 pam-fi 1.2 if(ncn.eq.0)then
196 pam-fi 1.8 if(debug.eq.1)print*,' WARNING - cnoise: ',
197 pam-fi 1.2 $ 'no strips for CN computation on VA1 ',j,
198 pam-fi 1.3 $ ', VIEW ',i,' >>> FAILED '
199 pam-fi 1.2 else
200 pam-fi 1.8 if(debug.eq.1)print*,' WARNING - cnoise: ',
201 pam-fi 1.3 $ 'less than ',NSTRIPMIN
202     $ ,' strips for CN computation on VA1 ',j,
203     $ ', VIEW ',i,' >>> FAILED '
204 pam-fi 1.2 endif
205 mocchiut 1.1 gulp=1
206 pam-fi 1.4 cnn(i,j) = 0
207 mocchiut 1.1 else
208 mocchiut 1.10 cn(i,j)=cn(i,j)/REAL(ncn) !<<<< computes common noise EM GCC4.7 CN IS REAL NOT DOUBLE
209     cnrms(i,j)= SQRT( cnrms(i,j)/REAL(ncn) - cn(i,j)**2 ) ! EM GCC4.7 CN IS REAL NOT DOUBLE
210 pam-fi 1.4 cnn(i,j) = ncn
211 pam-fi 1.2 gulp=0
212 pam-fi 1.7 c$$$ print*,'Event ',eventn(1)
213     c$$$ $ ,': cn(',i,',',j,')= ',cn(i,j),' ncn ',ncn
214    
215 pam-fi 1.8 if(debug.eq.1.and.ABS(cn(i,j)).gt.1000)
216 pam-fi 1.4 $ print*,'Event ',eventn(1)
217     $ ,': cn(',i,',',j,')= ',cn(i,j),' ncn ',ncn
218 mocchiut 1.1 endif
219    
220     return
221     end
222    
223    
224     *************************************************************************
225     *
226     * Subroutine cutcn.f!DA COMMENTARE!???
227     *
228     * excludes strips with particle signals and/or noisy strips from common
229     * noise calculation, marking their clstr(nviews,nva1_view,nstrips_va1)
230     * flag:
231     * clstr=0 ---> not to be used in CN computation
232     * clstr=1 ---> to be used in CN computation
233     *
234     * needs:
235     * - ./common_calib.f
236     *
237     * to be called inside ./cncomp.f
238     *
239     *************************************************************************
240    
241     subroutine cutcn(i,j) !(view, VA1)
242    
243     include 'commontracker.f'
244 pam-fi 1.4 include 'level1.f'
245 mocchiut 1.1 include 'common_reduction.f'
246     include 'calib.f'
247    
248    
249     integer skip !used to skip strips (see later...)
250    
251     integer kr, kl !position indexes to check signal affected
252     ! strips on right and left side of cluster
253     ! seed
254     integer ir, il !flags to exit loop on reaching VA1 extremes
255    
256 mocchiut 1.10 real(8) valuec !cluster seed signal EM GCC4.7
257 mocchiut 1.1 real cut,stripcut !cluster seed cut
258    
259 mocchiut 1.10 real(8) valuel, valuer !left and right strips signal EM GCC4.7
260 mocchiut 1.1 real stripcnincut !strip include cut
261    
262     skip = 0 !initializes skip
263    
264     do k=1,nstrips_va1 !loops on strips searching for cluster seeds
265    
266     if(k.le.skip) goto 20 !continues only if k strip has not been
267     ! checked yet
268    
269     clstr(i,j,k)=1 !reinitializes strip to be used in CN!???
270     ! computation, in order to be able to exclude
271     ! different strips at every CN computation loop
272    
273     c------------------------------------------------------------------------
274     c
275     c selects cut according to view
276     c
277     c------------------------------------------------------------------------
278     if(mod(i,2).eq.1) then !odd strip ---> Y view
279     valuec= - (DBLE(adc(i,j,k))-cn(i,j)-pedestal(i,j,k)) !negative signal
280     cut=clcuty !sets Y cut to find cluster seeds
281     else !even strip ---> X view
282     valuec= DBLE(adc(i,j,k))-cn(i,j)-pedestal(i,j,k) !positive signal
283     cut=clcutx !sets X cut to find cluster seeds
284     endif
285    
286    
287     c------------------------------------------------------------------------
288     c
289     c seeks clusters
290     c
291     c------------------------------------------------------------------------
292     stripcut=cut*sigma(i,j,k) !cluster seed cut
293    
294     c if(ABS(valuec).gt.stripcut) then !checks if signal exceeds threshold!???
295     if(valuec.gt.stripcut) then !checks if signal exceeds threshold
296    
297     c$$$ print*,'cut',i,j,k,valuec,stripcut,adc(i,j,k),cn(i,j)
298     c$$$ $ ,pedestal(i,j,k) !???
299    
300     clstr(i,j,k)=0 !if so, marks this strip as part of a
301     ! cluster
302    
303     c------------------------------------------------------------------------
304     c after finding a cluster seed, checks also adiacent strips, and marks
305     c the ones exceeding cnincut
306     c------------------------------------------------------------------------
307     kr=k !initializes position indexes to be equal to
308     kl=k ! cluster seed position
309    
310     ir=0 !initialize flags used to exit from
311     il=0 ! inclusion loop
312    
313     do while (il.eq.0.or.ir.eq.0) !shifts left and right from
314     ! cluster seed till it finds a strip below
315     ! the threshold, or till it reaches first or
316     ! last VA1 strip
317     kr=kr+1 !position index for strips on right side of
318     ! cluster seed
319     kl=kl-1 !and for left side
320    
321     c------------------------------------------------------------------------
322     c checks for last or first strip
323     c------------------------------------------------------------------------
324     if(kr.gt.nstrips_va1.and.ir.eq.0) then !when index goes
325     ir=1 ! beyond last VA1 strip, change ir flag in
326     ! order to "help" exiting from loop
327     skip=nstrips_va1+1 !sets skip beyond last strip: all
328     ! strips on the right have been included in
329     ! the cluster, so skips all next strips
330     ! (goto 20 condition is now always true)
331     endif
332    
333     if(kl.lt.1.and.il.eq.0) then !idem when index goes beyond
334     il=1 ! first strip
335     endif
336    
337     c P.S.: the "....and.i#.eq.0" term in above conditions is needed. In
338     c fact, even if I reach a under-cut strip on the right (so I get ir=1),
339     c the "do while loop" continues till such strip will be found on the
340     c left too.
341     c Thus kl and kr (! too) keep increasing, and it can happen kr gets
342     c greater than nstrips_va1 before kl reaches a under-cut strip. In this
343     c case it would pass this "if condition", so setting skip=nstrips_va1+1
344     c and skipping right strips never checked, if the "....and.i#.eq.0" term
345     c weren't the: instead, including this part it won't pass it
346     c because when I found reach the last VA1 strip on the right I set ir=1.
347     c (AAAAAAHHHHHHHHH!!!!!!!!!!!)
348    
349     c------------------------------------------------------------------------
350     c marks strips exceeding inclusion cut
351     c------------------------------------------------------------------------
352     c for right strips (kr index)
353     if(ir.eq.0) then !if last strip or last over-cut strip has
354     ! not been reached
355    
356     if(mod(i,2).eq.1) then !Y view
357     valuer= - (DBLE(adc(i,j,kr))-cn(i,j) !puts in valuer
358     $ -pedestal(i,j,kr)) ! right strip value
359     else !X view
360     valuer=DBLE(adc(i,j,kr))-cn(i,j)-pedestal(i,j,kr)
361     endif
362    
363     stripcnincut=cnincut*sigma(i,j,kr) !defines include cut
364     c if(ABS(valuer).gt.stripcnincut) then !marks right strip if it !???
365     if(valuer.gt.stripcnincut) then !marks right strip if it
366     clstr(i,j,kr)=0 !exceedes include cut
367     c$$$ print*,'inclcut_r',i,j,kr,valuer,stripcnincut
368     c$$$ $ ,adc(i,j,kr),cn(i,j),pedestal(i,j,kr) !???
369     else
370     ir=1 !otherwise cluster ends and ir flag =1
371     ! signals it
372     skip=kr !putting skip=kr, next k to be checked is
373     ! k=kr
374     endif
375    
376     endif
377    
378     c for left strips (kl index)
379     if(il.eq.0) then !if first strip or last over-cut strip has
380     ! not been reached
381    
382     if (mod(i,2).eq.1) then !Y view
383     valuel= - (DBLE(adc(i,j,kl))-cn(i,j) !puts in valuel
384     $ -pedestal(i,j,kl)) ! left strip value
385     else !X view
386     valuel=DBLE(adc(i,j,kl))-cn(i,j)-pedestal(i,j,kl)
387     endif
388    
389     stripcnincut=cnincut*sigma(i,j,kl) !defines include cut
390     c if(ABS(valuel).gt.stripcnincut) then !marks left strip if it!???
391     if(valuel.gt.stripcnincut) then !marks left strip if it
392     clstr(i,j,kl)=0 !exceedes include cut
393     c$$$ print*,'inclcut_l',i,j,kl,valuel,stripcnincut
394     c$$$ $ ,adc(i,j,kl),cn(i,j),pedestal(i,j,kl) !???
395     else
396     il=1 !otherwise cluster ends and il flag =1
397     ! signals it
398     endif
399    
400     endif
401    
402     enddo !ends lateral strips loop
403    
404     endif !ends cluster seed condition
405    
406     20 continue !comes here if next strip on the right has
407     ! already been included in a cluster
408    
409     enddo !ends principal strip loop
410    
411     return
412     end

  ViewVC Help
Powered by ViewVC 1.1.23