/[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.7 - (hide annotations) (download)
Thu May 24 16:45:48 2007 UTC (17 years, 6 months ago) by pam-fi
Branch: MAIN
Changes since 1.6: +3 -0 lines
several upgrades

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     if(debug)print*,' WARNING - cnoise: ',
175     $ 'no strips for CN computation on VA1 ',j,
176 pam-fi 1.3 $ ', VIEW ',i,' >>> FAILED '
177 pam-fi 1.2 else
178     if(debug)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.5 if(debug.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

  ViewVC Help
Powered by ViewVC 1.1.23