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