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