/[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.9 - (hide annotations) (download)
Fri Aug 31 14:56:52 2007 UTC (17 years, 3 months ago) by pam-fi
Branch: MAIN
CVS Tags: v5r00, v4r00, v9r00, v9r01, v6r01, v6r00
Changes since 1.8: +17 -0 lines
new variables added to TrkTrack + other changes

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

  ViewVC Help
Powered by ViewVC 1.1.23