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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.10 - (show 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
Error occurred while calculating annotation data.
Compilation warnings using GCC4.7 fixed

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 subroutine cncomp(i,j,errflag) !(view, VA1)
11
12 include 'commontracker.f'
13 include 'level1.f'
14 include 'common_reduction.f'
15 include 'calib.f'
16
17
18 integer errflag !error flag to mark no signal free VA1
19 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 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 integer newclstr !flag to warn about new found clusters to be
27 c ! excluded from common noise computation
28
29 c------------------------------------------------------------------------
30 c
31 c variables initialization
32 c
33 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 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 newclstr=1 !flag to warn about new found signal
58 c ! affected strips
59 c------------------------------------------------------------------------
60 c
61 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 c
65 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 nstr = nstr + strange(i,j,k) !uses only
75 if(mod(i,2).eq.1) then ! ---> Y view
76 signal(k) = - (DBLE(adc(i,j,k)) - pedestal(i,j,k)) !negative signal
77 else ! ---> X view
78 signal(k) = DBLE(adc(i,j,k)) - pedestal(i,j,k) !positive signal
79 endif
80 smean = smean + signal(k)*strange(i,j,k)
81 ssigma = ssigma + (signal(k)**2)*strange(i,j,k)
82 enddo
83
84 smean=smean/nstr !strips value distribution mean
85 ssigma=SQRT((ssigma/nstr)-smean**2) !strips value distribution sigma
86
87 cut=scut*ssigma !exclusion cut
88
89 nco=0
90 nbo=0
91 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 c print*,i,j,k,signal(k),smean
95 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
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 !???
110 if (countme.le.3) goto 666 !???
111
112 c------------------------------------------------------------------------
113 c
114 c common noise computation
115 c
116 c-----------------------------------------------------------------------
117 * loops on this VA1 till no new signal affected strips are found
118 do while(newclstr.eq.1)
119
120 newclstr=0 !to exit from loop if no new cluster is found
121
122 errflag=0
123 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
126 call cutcn(i,j) !(view, VA1) excludes clusters from cn computation
127
128 ncs=0 !initializes number of strips not excluded by cncut
129 do k=1,nstrips_va1 !loops on strips
130 * 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 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 enddo !ends do while
140
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 include 'level0.f'
171 include 'level1.f'
172 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 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
184 * tags strange, bad or signal-affected strips
185 iok = strange(i,j,k)*bad(i,j,k)*clstr(i,j,k)
186 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 ncn = ncn + iok !counts number of strips in cn computation
192 enddo
193
194 if(ncn.lt.NSTRIPMIN) then !no signal free strips on this VA1...
195 if(ncn.eq.0)then
196 if(debug.eq.1)print*,' WARNING - cnoise: ',
197 $ '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
206 cnn(i,j) = 0
207 else
208 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 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
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 include 'level1.f'
245 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 real(8) valuec !cluster seed signal EM GCC4.7
257 real cut,stripcut !cluster seed cut
258
259 real(8) valuel, valuer !left and right strips signal EM GCC4.7
260 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