/[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.6 - (hide annotations) (download)
Thu Oct 26 16:22:37 2006 UTC (18 years, 1 month ago) by pam-fi
Branch: MAIN
CVS Tags: v2r01, v3r04, v3r05, v3r06, v3r00, v3r03
Changes since 1.5: +8 -4 lines
fitting methods

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.5 if(debug.and.ABS(cn(i,j)).gt.1000)
191 pam-fi 1.4 $ print*,'Event ',eventn(1)
192     $ ,': cn(',i,',',j,')= ',cn(i,j),' ncn ',ncn
193 mocchiut 1.1 endif
194    
195     return
196     end
197    
198    
199     *************************************************************************
200     *
201     * Subroutine cutcn.f!DA COMMENTARE!???
202     *
203     * excludes strips with particle signals and/or noisy strips from common
204     * noise calculation, marking their clstr(nviews,nva1_view,nstrips_va1)
205     * flag:
206     * clstr=0 ---> not to be used in CN computation
207     * clstr=1 ---> to be used in CN computation
208     *
209     * needs:
210     * - ./common_calib.f
211     *
212     * to be called inside ./cncomp.f
213     *
214     *************************************************************************
215    
216     subroutine cutcn(i,j) !(view, VA1)
217    
218     include 'commontracker.f'
219 pam-fi 1.4 include 'level1.f'
220 mocchiut 1.1 include 'common_reduction.f'
221     include 'calib.f'
222    
223    
224     integer skip !used to skip strips (see later...)
225    
226     integer kr, kl !position indexes to check signal affected
227     ! strips on right and left side of cluster
228     ! seed
229     integer ir, il !flags to exit loop on reaching VA1 extremes
230    
231     real valuec !cluster seed signal
232     real cut,stripcut !cluster seed cut
233    
234     real valuel, valuer !left and right strips signal
235     real stripcnincut !strip include cut
236    
237     skip = 0 !initializes skip
238    
239     do k=1,nstrips_va1 !loops on strips searching for cluster seeds
240    
241     if(k.le.skip) goto 20 !continues only if k strip has not been
242     ! checked yet
243    
244     clstr(i,j,k)=1 !reinitializes strip to be used in CN!???
245     ! computation, in order to be able to exclude
246     ! different strips at every CN computation loop
247    
248     c------------------------------------------------------------------------
249     c
250     c selects cut according to view
251     c
252     c------------------------------------------------------------------------
253     if(mod(i,2).eq.1) then !odd strip ---> Y view
254     valuec= - (DBLE(adc(i,j,k))-cn(i,j)-pedestal(i,j,k)) !negative signal
255     cut=clcuty !sets Y cut to find cluster seeds
256     else !even strip ---> X view
257     valuec= DBLE(adc(i,j,k))-cn(i,j)-pedestal(i,j,k) !positive signal
258     cut=clcutx !sets X cut to find cluster seeds
259     endif
260    
261    
262     c------------------------------------------------------------------------
263     c
264     c seeks clusters
265     c
266     c------------------------------------------------------------------------
267     stripcut=cut*sigma(i,j,k) !cluster seed cut
268    
269     c if(ABS(valuec).gt.stripcut) then !checks if signal exceeds threshold!???
270     if(valuec.gt.stripcut) then !checks if signal exceeds threshold
271    
272     c$$$ print*,'cut',i,j,k,valuec,stripcut,adc(i,j,k),cn(i,j)
273     c$$$ $ ,pedestal(i,j,k) !???
274    
275     clstr(i,j,k)=0 !if so, marks this strip as part of a
276     ! cluster
277    
278     c------------------------------------------------------------------------
279     c after finding a cluster seed, checks also adiacent strips, and marks
280     c the ones exceeding cnincut
281     c------------------------------------------------------------------------
282     kr=k !initializes position indexes to be equal to
283     kl=k ! cluster seed position
284    
285     ir=0 !initialize flags used to exit from
286     il=0 ! inclusion loop
287    
288     do while (il.eq.0.or.ir.eq.0) !shifts left and right from
289     ! cluster seed till it finds a strip below
290     ! the threshold, or till it reaches first or
291     ! last VA1 strip
292     kr=kr+1 !position index for strips on right side of
293     ! cluster seed
294     kl=kl-1 !and for left side
295    
296     c------------------------------------------------------------------------
297     c checks for last or first strip
298     c------------------------------------------------------------------------
299     if(kr.gt.nstrips_va1.and.ir.eq.0) then !when index goes
300     ir=1 ! beyond last VA1 strip, change ir flag in
301     ! order to "help" exiting from loop
302     skip=nstrips_va1+1 !sets skip beyond last strip: all
303     ! strips on the right have been included in
304     ! the cluster, so skips all next strips
305     ! (goto 20 condition is now always true)
306     endif
307    
308     if(kl.lt.1.and.il.eq.0) then !idem when index goes beyond
309     il=1 ! first strip
310     endif
311    
312     c P.S.: the "....and.i#.eq.0" term in above conditions is needed. In
313     c fact, even if I reach a under-cut strip on the right (so I get ir=1),
314     c the "do while loop" continues till such strip will be found on the
315     c left too.
316     c Thus kl and kr (! too) keep increasing, and it can happen kr gets
317     c greater than nstrips_va1 before kl reaches a under-cut strip. In this
318     c case it would pass this "if condition", so setting skip=nstrips_va1+1
319     c and skipping right strips never checked, if the "....and.i#.eq.0" term
320     c weren't the: instead, including this part it won't pass it
321     c because when I found reach the last VA1 strip on the right I set ir=1.
322     c (AAAAAAHHHHHHHHH!!!!!!!!!!!)
323    
324     c------------------------------------------------------------------------
325     c marks strips exceeding inclusion cut
326     c------------------------------------------------------------------------
327     c for right strips (kr index)
328     if(ir.eq.0) then !if last strip or last over-cut strip has
329     ! not been reached
330    
331     if(mod(i,2).eq.1) then !Y view
332     valuer= - (DBLE(adc(i,j,kr))-cn(i,j) !puts in valuer
333     $ -pedestal(i,j,kr)) ! right strip value
334     else !X view
335     valuer=DBLE(adc(i,j,kr))-cn(i,j)-pedestal(i,j,kr)
336     endif
337    
338     stripcnincut=cnincut*sigma(i,j,kr) !defines include cut
339     c if(ABS(valuer).gt.stripcnincut) then !marks right strip if it !???
340     if(valuer.gt.stripcnincut) then !marks right strip if it
341     clstr(i,j,kr)=0 !exceedes include cut
342     c$$$ print*,'inclcut_r',i,j,kr,valuer,stripcnincut
343     c$$$ $ ,adc(i,j,kr),cn(i,j),pedestal(i,j,kr) !???
344     else
345     ir=1 !otherwise cluster ends and ir flag =1
346     ! signals it
347     skip=kr !putting skip=kr, next k to be checked is
348     ! k=kr
349     endif
350    
351     endif
352    
353     c for left strips (kl index)
354     if(il.eq.0) then !if first strip or last over-cut strip has
355     ! not been reached
356    
357     if (mod(i,2).eq.1) then !Y view
358     valuel= - (DBLE(adc(i,j,kl))-cn(i,j) !puts in valuel
359     $ -pedestal(i,j,kl)) ! left strip value
360     else !X view
361     valuel=DBLE(adc(i,j,kl))-cn(i,j)-pedestal(i,j,kl)
362     endif
363    
364     stripcnincut=cnincut*sigma(i,j,kl) !defines include cut
365     c if(ABS(valuel).gt.stripcnincut) then !marks left strip if it!???
366     if(valuel.gt.stripcnincut) then !marks left strip if it
367     clstr(i,j,kl)=0 !exceedes include cut
368     c$$$ print*,'inclcut_l',i,j,kl,valuel,stripcnincut
369     c$$$ $ ,adc(i,j,kl),cn(i,j),pedestal(i,j,kl) !???
370     else
371     il=1 !otherwise cluster ends and il flag =1
372     ! signals it
373     endif
374    
375     endif
376    
377     enddo !ends lateral strips loop
378    
379     endif !ends cluster seed condition
380    
381     20 continue !comes here if next strip on the right has
382     ! already been included in a cluster
383    
384     enddo !ends principal strip loop
385    
386     return
387     end

  ViewVC Help
Powered by ViewVC 1.1.23