/[PAMELA software]/DarthVader/TrackerLevel2/src/F77/reductionflight.f
ViewVC logotype

Annotation of /DarthVader/TrackerLevel2/src/F77/reductionflight.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.26 - (hide annotations) (download)
Tue Nov 25 14:41:38 2008 UTC (16 years ago) by pam-fi
Branch: MAIN
Changes since 1.25: +31 -9 lines
fixed small bug in cluster-finding.

1 mocchiut 1.1 *************************************************************************
2     *
3     * Program reductionflight.f
4     *
5     * - reads readraw.f output files: LEVEL0 ntuple, and ped, sig and bad histograms
6     * - decodes raw data (DATATRACKER) using DSP ped, sig and bad values
7     * - looks for clusters information using ped, sig and bad values from
8     * DSP histograms
9     * - fills LEVEL1 ntuple
10     *
11     *************************************************************************
12    
13 pam-fi 1.2 subroutine reductionflight(ierror)
14 mocchiut 1.1
15     include 'commontracker.f'
16     include 'level0.f'
17     include 'level1.f'
18     include 'common_reduction.f'
19     include 'calib.f'
20    
21 pam-fi 1.6 data eventn_old/nviews*0/
22    
23 pam-fi 1.2 integer ierror
24     ierror = 0
25 mocchiut 1.1
26 pam-fi 1.17 c$$$ debug = .true.
27     c$$$ verbose = .true.
28     c$$$ warning = .true.
29    
30 pam-fi 1.22 c$$$ print*,debug,verbose,warning
31     c$$$ debug=1
32     c$$$ verbose=1
33     c$$$ warning=1
34    
35 pam-fi 1.19 * //////////////////////////
36     * initialize some parameters
37     * //////////////////////////
38    
39 mocchiut 1.1 call init_level1
40    
41 pam-fi 1.20 c debug=.true.
42    
43 pam-fi 1.22 if(debug.eq.1)print*,'-- check LEVEL0 status'
44 pam-fi 1.17
45 pam-fi 1.20 ievco=-1
46     mismatch=0
47 pam-fi 1.6 c good1 = good0
48     c--------------------------------------------------
49     c check the LEVEL0 event status for missing
50     c sections or DSP alarms
51     c ==> set the variable GOOD1(12)
52     c--------------------------------------------------
53     do iv=1,nviews
54     if(DSPnumber(iv).gt.0.and.DSPnumber(iv).le.12)then
55     c ------------------------
56     c GOOD
57     c ------------------------
58 pam-fi 1.9 GOOD1(DSPnumber(iv))=0 !OK
59 pam-fi 1.6 c ------------------------
60     c CRC error
61     c ------------------------
62     if(crc(iv).eq.1) then
63 pam-fi 1.20 c GOOD1(DSPnumber(iv)) = 2
64     c GOOD1(DSPnumber(iv)) = GOOD1(DSPnumber(iv)) + 2**1
65     GOOD1(DSPnumber(iv)) = ior(GOOD1(DSPnumber(iv)),2**1)
66     102 format(' * WARNING * Event ',i7,' view',i3
67     $ ,' CRC error')
68 pam-fi 1.22 if(debug.eq.1)write(*,102)eventn(1),DSPnumber(iv)
69 pam-fi 1.20 c goto 18 !next view
70 pam-fi 1.6 endif
71     c ------------------------
72     c online-software alarm
73     c ------------------------
74     if(
75     $ fl1(iv).ne.0.or.
76     $ fl2(iv).ne.0.or.
77     $ fl3(iv).ne.0.or.
78     $ fl4(iv).ne.0.or.
79     $ fl5(iv).ne.0.or.
80     $ fl6(iv).ne.0.or.
81     $ fc(iv).ne.0.or.
82     $ DATAlength(iv).eq.0.or.
83     $ .false.)then
84 pam-fi 1.20 c GOOD1(DSPnumber(iv))=3
85     c GOOD1(DSPnumber(iv)) = GOOD1(DSPnumber(iv)) + 2**2
86     GOOD1(DSPnumber(iv)) = ior(GOOD1(DSPnumber(iv)),2**2)
87     103 format(' * WARNING * Event ',i7,' view',i3
88     $ ,' software alarm')
89 pam-fi 1.22 if(debug.eq.1)write(*,103)eventn(1),DSPnumber(iv)
90 pam-fi 1.20 c goto 18
91 pam-fi 1.6 endif
92     c ------------------------
93     c DSP-counter jump
94     c ------------------------
95 pam-fi 1.20 c commentato perche` non e` un controllo significativo nel caso in cui
96     c la subroutine venga chiamata per riprocessare l'evento
97     c sostituito con un check dei contatori dei vari dsp
98     c$$$ if(
99     c$$$ $ eventn_old(iv).ne.0.and. !first event in this file
100     c$$$ $ eventn(iv).ne.1.and. !first event in run
101     c$$$ $ good_old(DSPnumber(iv)).ne.0.and. !previous event corrupted
102     c$$$ $ .true.)then
103     c$$$
104     c$$$ if(eventn(iv).ne.(eventn_old(iv)+1))then
105     c$$$c GOOD1(DSPnumber(iv))=4
106     c$$$c GOOD1(DSPnumber(iv)) = GOOD1(DSPnumber(iv)) + 2**3
107     c$$$ GOOD1(DSPnumber(iv)) = ior(GOOD1(DSPnumber(iv)),2**3)
108     c$$$ 104 format(' * WARNING * Event ',i7,' view',i3
109     c$$$ $ ,' counter jump ',i10,i10)
110 pam-fi 1.22 c$$$ if(debug.eq.1)write(*,104)eventn(1),DSPnumber(iv)
111 pam-fi 1.20 c$$$ $ ,eventn_old(iv),eventn(iv))
112     c$$$ goto 18
113     c$$$ endif
114     c$$$
115     c$$$ endif
116     c ------------------------
117     c 18 continue
118     c ------------------------
119     c DSP-counter
120     c ------------------------
121     if( DSPnumber(iv).ne.0.and.GOOD1(DSPnumber(iv)).ne.1)then
122     if(iv.ne.1.and.ievco.ne.-1)then
123     if( eventn(iv).ne.ievco )then
124     mismatch=1
125     endif
126 pam-fi 1.6 endif
127 pam-fi 1.20 ievco = eventn(iv)
128 pam-fi 1.6 endif
129     endif
130     enddo
131    
132 pam-fi 1.20 c print*,'*** ',(eventn(iv),iv=1,12)
133    
134 pam-fi 1.22 if(mismatch.eq.1.and.debug.eq.1)
135 pam-fi 1.20 $ print*,' * WARNING * DSP counter mismatch: '
136     $ ,(eventn(iv),iv=1,12)
137    
138 pam-fi 1.6 ngood = 0
139     do iv = 1,nviews
140 pam-fi 1.20
141     if(mismatch.eq.1.and.GOOD1(iv).ne.1)
142     $ GOOD1(iv)=ior(GOOD1(iv),2**3)
143    
144 pam-fi 1.6 eventn_old(iv) = eventn(iv)
145     good_old(iv) = good1(iv)
146     ngood = ngood + good1(iv)
147 pam-fi 1.20
148 pam-fi 1.6 enddo
149 pam-fi 1.20 c$$$ if(verbose.and.ngood.ne.0)print*,'* WARNING * Event ',eventn(1)
150     c$$$ $ ,':LEVEL0 event status: '
151     c$$$ $ ,(good1(i),i=1,nviews)
152 mocchiut 1.1 c--------------------------------------------------
153     c read the variable DATATRACKER from LEVEL0
154 pam-fi 1.6 c and fill the variable ADC (invertin view 11)
155 mocchiut 1.1 c--------------------------------------------------
156 pam-fi 1.17
157 pam-fi 1.22 if(debug.eq.1)print*,'-- fill ADC vectors'
158 pam-fi 1.17
159 mocchiut 1.1 call filladc(iflag)
160     if(iflag.ne.0)then
161 pam-fi 1.6 ierror = 220
162 mocchiut 1.1 endif
163    
164     c--------------------------------------------------
165     c computes common noise for each VA1
166 pam-fi 1.10 c (excluding strips with signal,
167 mocchiut 1.1 c tagged with the flag CLSTR)
168     c--------------------------------------------------
169 pam-fi 1.22 if(debug.eq.1)print*,'-- compute CN'
170 pam-fi 1.17
171 mocchiut 1.1 do iv=1,nviews
172 pam-fi 1.25
173     call evaluatecn(iv)
174     c$$$ ima=0
175     c$$$ do ik=1,nva1_view
176     c$$$ cn(iv,ik) = 0
177     c$$$ cnrms(iv,ik) = 0
178     c$$$ cnn(iv,ik) = -1
179     c$$$ iflag = 0
180     c$$$ mask_vk_ev(iv,ik) = 1
181     c$$$ call stripmask(iv,ik) !compute mask(i,j,k), combining VA1-masks
182     c$$$* --------------------------------------
183     c$$$* if chip is not masked ---> evaluate CN
184     c$$$* --------------------------------------
185     c$$$ if( mask(iv,ik,1).eq.1 ) then !!!NBNB mask per la striscia 1 !!!!!!!!
186     c$$$ call cncomp(iv,ik,iflag)
187     c$$$ if(iflag.ne.0)then
188     c$$$ ima=ima+1
189     c$$$ mask_vk_ev(iv,ik)=0
190     c$$$ ierror = 220
191     c$$$ endif
192     c$$$ call stripmask(iv,ik) !compute mask(i,j,k), combining VA1-masks
193     c$$$ endif
194     c$$$ enddo
195     c$$$ 100 format(' * WARNING * Event ',i7,' view',i3,': VK MASK ',24i1)
196     c$$$ if(ima.ne.0.and.verbose.eq.1)write(*,100)eventn(1),iv
197     c$$$ $ ,(mask_vk_ev(iv,ik),ik=1,nva1_view)
198    
199 mocchiut 1.1 enddo
200    
201 pam-fi 1.13 cc call stripmask !compute mask(i,j,k), combining mask_vk_ev and mask_vk
202 pam-fi 1.10
203 mocchiut 1.1 c---------------------------------------------
204     c loops on views, VA1 and strips,
205     c and computes strips signals using
206     c badstrip, pedestals, and
207     c sigma informations from histograms
208     c---------------------------------------------
209     ind=1 !clsignal array index
210 pam-fi 1.5
211 pam-fi 1.22 if(debug.eq.1)print*,'-- search clusters'
212 mocchiut 1.1 do iv=1,nviews !loop on views
213 pam-fi 1.25 c$$$ do is=1,nstrips_view !loop on strips (1)
214     c$$$ if(mod(iv,2).eq.1) then
215     c$$$C=== > Y view
216     c$$$c print*,iv,nvk(is),nst(is),adc(iv,nvk(is),nst(is))
217     c$$$c $ ,cn(iv,nvk(is))
218     c$$$c $ ,pedestal(iv,nvk(is),nst(is))
219     c$$$ value(is)= -(DBLE(adc(iv,nvk(is),nst(is)))
220     c$$$ $ -cn(iv,nvk(is))-pedestal(iv,nvk(is),nst(is)))
221     c$$$ $ *mask(iv,nvk(is),nst(is))
222     c$$$ clseedcut(is)=clcuty*sigma(iv,nvk(is),nst(is))
223     c$$$ $ *mask(iv,nvk(is),nst(is))
224     c$$$ clinclcut(is)=incuty*sigma(iv,nvk(is),nst(is))
225     c$$$ $ *mask(iv,nvk(is),nst(is))
226     c$$$ sat(is)=0
227     c$$$ if( adc(iv,nvk(is),nst(is)).lt.adc_saty )sat(is)=1
228     c$$$ else
229     c$$$C=== > X view
230     c$$$ value(is)= (DBLE(adc(iv,nvk(is),nst(is)))
231     c$$$ $ -cn(iv,nvk(is))-pedestal(iv,nvk(is),nst(is)))
232     c$$$ $ *mask(iv,nvk(is),nst(is))
233     c$$$ clseedcut(is)=clcutx*sigma(iv,nvk(is),nst(is))
234     c$$$ $ *mask(iv,nvk(is),nst(is))
235     c$$$ clinclcut(is)=incutx*sigma(iv,nvk(is),nst(is))
236     c$$$ $ *mask(iv,nvk(is),nst(is))
237     c$$$ sat(is)=0
238     c$$$ if( adc(iv,nvk(is),nst(is)).gt.adc_satx )sat(is)=1
239     c$$$ endif
240     c$$$ enddo !end loop on strips (1)
241     call subtractped(iv)
242     call searchcluster(iv)
243    
244     if(.not.flag_shower)then
245     call savecluster(iv)
246 pam-fi 1.26 c$$$ if(debug.eq.1)print*,' view ',iv,' #clusters ', nclstr_view
247 pam-fi 1.25 else
248     fshower(iv) = 1
249     c GOOD1(DSPnumber(iv)) = 11 !AHAHAHAHA ORRORE!!
250     c GOOD1(iv) = 11
251     c GOOD1(iv) = GOOD1(iv) + 2**5
252     GOOD1(iv) = ior(GOOD1(iv),2**5)
253     101 format(' * WARNING * Event ',i7,' view',i3
254     $ ,' #clusters > ',i5,' --> MASKED')
255     if(verbose.eq.1)write(*,101)eventn(1),iv,nclstrmax_view
256     endif
257 mocchiut 1.1 enddo ! end loop on views
258     do iv=1,nviews
259 pam-fi 1.25 do ik=1,nva1_view
260     cnev(iv,ik) = cn(iv,ik) !assigns computed CN to ntuple variables
261     cnrmsev(iv,ik) = cnrms(iv,ik) !assigns computed CN to ntuple variables
262     cnnev(iv,ik) = cnn(iv,ik) !assigns computed CN to ntuple variables
263     enddo
264 mocchiut 1.1 enddo
265     C---------------------------------------------
266     C come here if GOOD1=0
267     C or the event has too many clusters
268     C---------------------------------------------
269     200 continue
270 pam-fi 1.6
271     ngood = 0
272     do iv = 1,nviews
273     ngood = ngood + good1(iv)
274     enddo
275 pam-fi 1.22 if(verbose.eq.1.and.ngood.ne.0)
276     $ print*,'* WARNING * Event ',eventn(1)
277 pam-fi 1.17 $ ,':LEVEL1 event status: '
278     $ ,(good1(i),i=1,nviews)
279 mocchiut 1.1 c------------------------------------------------------------------------
280 pam-fi 1.2 c
281 mocchiut 1.1 c closes files and exits
282 pam-fi 1.2 c
283 mocchiut 1.1 c------------------------------------------------------------------------
284 pam-fi 1.2 RETURN
285     END
286 mocchiut 1.1
287     ***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***
288     *
289     *
290     *
291     *
292     *
293     *
294     *
295     *
296     *
297     ***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***
298    
299    
300     subroutine init_level1
301    
302     include 'commontracker.f'
303     include 'level1.f'
304     include 'level0.f'
305    
306 pam-fi 1.6 c good1 = 0
307     do iv=1,12
308     good1(iv) = 1 !missing packet
309     enddo
310 pam-fi 1.5 nclstr1 = 0
311     totCLlength = 0
312 mocchiut 1.1 do ic=1,nclstrmax
313 pam-fi 1.5 view(ic) = 0
314     ladder(ic) = 0
315     indstart(ic) = 0
316     indmax(ic) = 0
317     maxs(ic) = 0
318     mult(ic) = 0
319 pam-fi 1.16 sgnl(ic) = 0
320 pam-fi 1.15 whichtrack(ic) = 0 !assigned @ level2
321 pam-fi 1.5
322 mocchiut 1.1 enddo
323     do id=1,maxlength !???
324 pam-fi 1.5 clsignal(id) = 0.
325     clsigma(id) = 0.
326     cladc(id) = 0.
327     clbad(id) = 0.
328 mocchiut 1.1 enddo
329     do iv=1,nviews
330     c crc1(iv)=0
331     do ik=1,nva1_view
332 pam-fi 1.5 cnev(iv,ik) = 0
333     cnnev(iv,ik) = 0
334 mocchiut 1.1 enddo
335 pam-fi 1.5 fshower(iv) = 0
336 mocchiut 1.1 enddo
337    
338     return
339     end
340 pam-fi 1.10
341 mocchiut 1.1 *---***---***---***---***---***---***---***---***
342     *
343     *
344     *
345     *
346     *
347     *---***---***---***---***---***---***---***---***
348    
349 pam-fi 1.25 subroutine searchcluster(iv)
350 mocchiut 1.1
351     include 'commontracker.f'
352     include 'level0.f'
353     include 'level1.f'
354     include 'calib.f'
355    
356 pam-fi 1.5 include 'common_reduction.f'
357 mocchiut 1.1
358    
359     c local variables
360     integer rmax,lmax !estremi del cluster
361 pam-fi 1.10 integer rstop,lstop
362     integer first,last
363     integer fsat,lsat
364 mocchiut 1.1
365     external nst
366    
367 pam-fi 1.10 iseed=-999 !cluster seed index initialization
368 mocchiut 1.1
369 pam-fi 1.10 inext=-999 !index where to start new cluster search
370 mocchiut 1.1
371 pam-fi 1.10 flag_shower = .false.
372 pam-fi 1.5 nclstr_view=0
373    
374 pam-fi 1.10 do jl=1,nladders_view !1..3 !loops on ladders
375    
376     first = 1+nstrips_ladder*(jl-1) !1,1025,2049
377     last = nstrips_ladder*jl !1024,2048,3072
378    
379     * X views have 1018 strips instead of 1024
380 mocchiut 1.1 if(mod(iv,2).eq.0) then
381     first=first+3
382     last=last-3
383     endif
384 pam-fi 1.6
385 mocchiut 1.1 do is=first,last !loop on strips in each ladder
386 pam-fi 1.6
387 pam-fi 1.10 c---------------------------------------------
388     c new-cluster search starts at index inext
389     c---------------------------------------------
390     if(is.lt.inext) goto 220 ! next strip
391 pam-fi 1.5
392 mocchiut 1.1 if(value(is).gt.clseedcut(is)) then
393     c-----------------------------------------
394     c possible SEED...
395     c-----------------------------------------
396 pam-fi 1.26 c$$$ if(debug.eq.1)print*,'|||| ',value(is),' @',is
397     c$$$ $ ,' cut ',clseedcut(is)
398    
399 pam-fi 1.10 itemp = is
400     fsat = 0 ! first saturated strip
401     lsat = 0 ! last saturated strip
402 mocchiut 1.1 if(itemp.eq.last) goto 230 !estremo...
403 pam-fi 1.10 c ------------------------
404     c search for first maximum
405     c ------------------------
406 pam-fi 1.6 do while(
407     $ value(itemp).le.value(itemp+1)
408     $ .and.value(itemp+1).gt.clseedcut(itemp+1))
409 pam-fi 1.10 itemp = itemp+1
410 pam-fi 1.26 c$$$ if(debug.eq.1)print*,'|||| ',value(itemp),' @',is
411     c$$$ $ ,' cut ',clseedcut(itemp)
412 pam-fi 1.10 if(itemp.eq.last) goto 230 !stops if reaches last strip
413     if(sat(itemp).eq.1) goto 230 !stop if reaches a saturated strip
414 mocchiut 1.1 enddo ! of the ladder
415     230 continue
416 pam-fi 1.10 c -----------------------------
417     c check if strips are saturated
418     c -----------------------------
419     if( sat(itemp).eq.1 )then
420     fsat = itemp
421     lsat = itemp
422     if(itemp.eq.last) goto 231 !estremo...
423 pam-fi 1.26 do while(
424     $ sat(itemp+1).eq.1 .and.
425     $ value(itemp+1).gt.clseedcut(itemp+1) .and.
426     $ .true.)
427 pam-fi 1.10 itemp = itemp+1
428     lsat = itemp
429     if(itemp.eq.last) goto 231 !stops if reaches last strip
430     enddo
431     endif
432     231 continue
433     c---------------------------------------------------------------------------
434 mocchiut 1.1 c fownd SEED!!!
435 pam-fi 1.10 c (if there are saturated strips, the cluster is centered in the middle)
436     c---------------------------------------------------------------------------
437     if(fsat.eq.0.and.lsat.eq.0)then
438     iseed = itemp ! <<< SEED
439     else
440     iseed = int((lsat+fsat)/2) ! <<< SEED
441 pam-fi 1.26 if(debug.eq.1)
442     $ print*,'saturated strips (first,last) ',fsat,lsat
443 pam-fi 1.10 c$$$ print*,'--> ',(value(ii),ii=fsat,lsat)
444     endif
445     c---------------------------------------------------------------
446 mocchiut 1.1 c after finding a cluster seed, checks also adjacent strips,
447 pam-fi 1.10 C and tags the ones exceeding clinclcut
448     c---------------------------------------------------------------
449 pam-fi 1.26
450     if(debug.eq.1)print*,'SEED ',value(iseed),' @',iseed
451     $ ,' cut ',clseedcut(iseed)
452    
453 mocchiut 1.1 ir=iseed !indici destro
454     il=iseed ! e sinistro
455    
456     rmax=ir !estremo destro del cluster
457     lmax=il ! e sinistro
458    
459     rstop=0 !initialize flags used to exit from
460     lstop=0 ! inclusion loop
461    
462     do while(lstop.eq.0.or.rstop.eq.0) !shifts left and right from
463 pam-fi 1.10
464    
465     ir=ir+1 !index for right side
466     il=il-1 !index for left side
467 mocchiut 1.1 c------------------------------------------------------------------------
468     c checks for last or first strip of the ladder
469     c------------------------------------------------------------------------
470 pam-fi 1.10 if( ir.gt.last ) rstop = 1
471     if( il.lt.first ) lstop = 1
472 mocchiut 1.1
473     c------------------------------------------------------------------------
474 pam-fi 1.10 c add strips exceeding inclusion cut
475 mocchiut 1.1 c------------------------------------------------------------------------
476 pam-fi 1.10 if( (rmax-lmax+1).ge.nclstrp )goto 210 !exits inclusion loop
477    
478 pam-fi 1.21 if(rstop.eq.0) then !if right cluster border has not been reached
479 pam-fi 1.10 if(value(ir).gt.clinclcut(ir)) then
480     rmax=ir !include a strip on the right
481 mocchiut 1.1 else
482 pam-fi 1.10 rstop=1 !cluster right end
483     endif
484 mocchiut 1.1 endif
485 pam-fi 1.10
486     if( (rmax-lmax+1).ge.nclstrp )goto 210 !exits inclusion loop
487    
488 pam-fi 1.21 if(lstop.eq.0) then !if left cluster border has not been reached
489 mocchiut 1.1 if(value(il).gt.clinclcut(il)) then
490 pam-fi 1.10 lmax=il !include a strip on the left
491 mocchiut 1.1 else
492 pam-fi 1.10 lstop=1 !cluster left end
493 mocchiut 1.1 endif
494     endif
495    
496 pam-fi 1.21 c if( (rmax-lmax+1).ge.nclstrp )goto 210 !exits inclusion loop
497    
498 mocchiut 1.1 enddo !ends strip inclusion loop
499 pam-fi 1.10 goto 211
500 mocchiut 1.1 210 continue !jumps here if more than nclstrp have been included
501 pam-fi 1.10 c print*,'>>> nclstrp! '
502     211 continue
503     c-----------------------------------------
504     c end of inclusion loop!
505     c-----------------------------------------
506    
507     c------------------------------------------------------------------------
508     c adjust the cluster in order to have at least a strip around the seed
509     c------------------------------------------------------------------------
510     if(iseed.eq.lmax.and.lmax.ne.first)then
511     lmax = lmax-1
512     if( (rmax-lmax+1).gt.nclstrp )rmax=rmax-1
513     endif
514     if(iseed.eq.rmax.and.rmax.ne.last )then
515     rmax = rmax+1
516     if( (rmax-lmax+1).gt.nclstrp )lmax=lmax+1
517 mocchiut 1.1 endif
518 pam-fi 1.21 c-------------------------------------------------------------------------------
519     c adjust the cluster in order to have at least ANOTHER strip around the seed
520     c-------------------------------------------------------------------------------
521 pam-fi 1.24 if(iseed-1.eq.lmax.and.lmax.ne.first)then
522     lmax = lmax-1
523     if( (rmax-lmax+1).gt.nclstrp )rmax=rmax-1
524     endif
525     if(iseed+1.eq.rmax.and.rmax.ne.last )then
526     rmax = rmax+1
527     if( (rmax-lmax+1).gt.nclstrp )lmax=lmax+1
528     endif
529 pam-fi 1.21 c---------------------------------------------------
530     c now we have 5 stored-strips around the maximum
531     c---------------------------------------------------
532 pam-fi 1.10
533 mocchiut 1.1 c------------------------------------------------------------------------
534 pam-fi 1.10 c adjust the cluster in order to store a minimum number of strips
535     c------------------------------------------------------------------------
536     do while( (rmax-lmax+1).lt.nclstrpmin )
537    
538     vl = -99999
539     vr = -99999
540     if(lmax-1.ge.first) vl = value(lmax-1)
541     if(rmax+1.le.last ) vr = value(rmax+1)
542     if(vr.ge.vl) then
543     rmax = rmax+1
544     else
545     lmax = lmax-1
546 mocchiut 1.1 endif
547 pam-fi 1.10
548     enddo
549 mocchiut 1.1
550     c--------------------------------------------------------
551 pam-fi 1.10 c store cluster info
552 mocchiut 1.1 c--------------------------------------------------------
553 pam-fi 1.5 nclstr_view = nclstr_view + 1 !cluster number
554 pam-fi 1.10
555 pam-fi 1.5 if(nclstr_view.gt.nclstrmax_view) then !too many clusters for the view:
556 pam-fi 1.10 c$$$ if(verbose) print*,'Event ',eventn(1),
557     c$$$ $ ': more than ',nclstrmax_view
558     c$$$ $ ,' clusters on view ',iv
559 mocchiut 1.1 flag_shower = .true.
560     goto 2000
561     endif
562 pam-fi 1.5
563 pam-fi 1.10 ladder_view(nclstr_view) = nld(iseed,iv)
564     maxs_view(nclstr_view) = iseed
565 pam-fi 1.5 rmax_view(nclstr_view) = rmax
566     lmax_view(nclstr_view) = lmax
567 pam-fi 1.23 c mult_view(nclstr_view) = rmax-lmax+1
568     mult_view(nclstr_view) = 0
569     do ii=lmax,rmax
570     if(value(ii).gt.clinclcut(ii))
571     $ mult_view(nclstr_view) = mult_view(nclstr_view)+1
572     enddo
573    
574 pam-fi 1.26 c$$$ print*,(value(ii),ii=lmax,rmax)
575     c$$$ print*,(clinclcut(ii),ii=lmax,rmax)
576     c$$$ print*,(clseedcut(ii),ii=lmax,rmax)
577 pam-fi 1.5
578 pam-fi 1.10 c$$$ if(rmax-lmax+1.gt.25)
579     c$$$ $ print*,'view ',iv
580     c$$$ $ ,' cl ',nclstr_view,' mult ',rmax-lmax+1
581     c------------------------------------------------------------------------
582 pam-fi 1.11 c search for a double peak inside the cluster
583 pam-fi 1.10 c------------------------------------------------------------------------
584     inext = rmax+1 !<< index where to start new-cluster search
585    
586     vmax = 0
587     vmin = value(iseed)
588     imax = iseed
589     imin = iseed
590     do iss = max(iseed+1,lsat+1),rmax
591     if( value(iss).lt.vmin )then
592     if( imax.ne.iseed )goto 221 !found dowble peek
593     imin = iss
594     vmin = value(iss)
595     else
596     delta = value(iss) - value(imin)
597     cut = sqrt(clinclcut(iss)**2 + clinclcut(imin)**2)
598     if(
599     $ delta.gt.cut .and.
600     $ value(iss).gt.clseedcut(iss).and.
601     $ .true.)then
602     if( value(iss).gt.vmax )then
603     imax = iss
604     vmax = value(iss)
605     else
606     goto 221 !found dowble peek
607     endif
608     endif
609     endif
610     enddo
611     221 continue
612    
613     if(imax.gt.iseed)then
614     inext = imax !<< index where to start new-cluster search
615     c$$$ print*,'--- double peek ---'
616     c$$$ print*,(value(ii),ii=lmax,rmax)
617     c$$$ print*,'seed ',iseed,' imin ',imin,' imax ',imax
618     endif
619 mocchiut 1.1 c--------------------------------------------------------
620 pam-fi 1.2 c
621 mocchiut 1.1 c--------------------------------------------------------
622     endif !end possible seed conditio
623     220 continue !jumps here to skip strips left of last seed
624    
625     enddo ! end loop on strips
626     enddo !end loop on ladders
627     2000 continue
628     return
629     end
630    
631    
632     *---***---***---***---***---***---***---***---***
633     *
634     *
635     *
636     *
637     *
638     *---***---***---***---***---***---***---***---***
639    
640 pam-fi 1.25 subroutine savecluster(iv)
641 pam-fi 1.5 *
642     * (080/2006 Elena Vannuccini)
643     * Save the clusters view by view
644    
645     include 'commontracker.f'
646     include 'level1.f'
647     include 'calib.f'
648     include 'common_reduction.f'
649    
650     integer CLlength !lunghezza in strip del cluster
651    
652     do ic=1,nclstr_view
653    
654     nclstr1 = nclstr1+1
655     view(nclstr1) = iv
656     ladder(nclstr1) = ladder_view(ic)
657     maxs(nclstr1) = maxs_view(ic)
658     mult(nclstr1) = mult_view(ic)
659    
660     c posizione dell'inizio del cluster nell' array clsignal
661     indstart(nclstr1) = ind
662     c posizione del cluster seed nell'array clsignal
663     indmax(nclstr1) = indstart(nclstr1)
664     $ +( maxs_view(ic) - lmax_view(ic) )
665    
666     CLlength = rmax_view(ic) - lmax_view(ic) + 1 !numero di strip salvate
667     totCLlength = totCLlength + CLlength
668 pam-fi 1.16 sgnl(nclstr1) = 0
669 pam-fi 1.5 do j=lmax_view(ic),rmax_view(ic) !stores sequentially cluter strip values in
670    
671     clsignal(ind) = value(j) ! clsignal array
672 pam-fi 1.20 c$$$ print*,ind,clsignal(ind)
673 pam-fi 1.5 ivk=nvk(j)
674     ist=nst(j)
675    
676     clsigma(ind) = sigma(iv,ivk,ist)
677     cladc(ind) = adc(iv,ivk,ist)
678     clbad(ind) = bad(iv,ivk,ist)
679     c clped(ind) = pedestal(iv,ivk,ist)
680    
681     ind=ind+1
682     c if(value(j).gt.0)
683     if(value(j).gt.clinclcut(j))
684 pam-fi 1.16 $ sgnl(nclstr1) = sgnl(nclstr1) + value(j) !cluster charge
685 pam-fi 1.5 enddo
686    
687 pam-fi 1.26 if(debug.eq.1)then
688     print*,'view ',iv,' -- '
689     $ ,' n.cl: ',nclstr1
690     $ ,' maxs: ',maxs(nclstr1)
691     $ ,' mult: ',mult(nclstr1)
692     $ ,' sign: ',sgnl(nclstr1)
693     print*,'----------------------'
694     endif
695 pam-fi 1.5 enddo
696    
697     return
698     end
699     *---***---***---***---***---***---***---***---***
700     *
701     *
702     *
703     *
704     *
705     *---***---***---***---***---***---***---***---***
706    
707 pam-fi 1.25 subroutine evaluatecn(iv)
708    
709     include 'commontracker.f'
710     include 'level0.f'
711     include 'level1.f'
712     include 'common_reduction.f'
713     include 'calib.f'
714    
715     ima=0
716     do ik=1,nva1_view
717     cn(iv,ik) = 0
718     cnrms(iv,ik) = 0
719     cnn(iv,ik) = -1
720     iflag = 0
721     mask_vk_ev(iv,ik) = 1
722     call stripmask(iv,ik) !compute mask(i,j,k), combining VA1-masks
723     * --------------------------------------
724     * if chip is not masked ---> evaluate CN
725     * --------------------------------------
726     if( mask(iv,ik,1).eq.1 ) then !!!NBNB mask per la striscia 1 !!!!!!!!
727     call cncomp(iv,ik,iflag)
728     if(iflag.ne.0)then
729     ima=ima+1
730     mask_vk_ev(iv,ik)=0
731     ierror = 220
732     endif
733     call stripmask(iv,ik) !compute mask(i,j,k), combining VA1-masks
734     endif
735     enddo
736     100 format(' * WARNING * Event ',i7,' view',i3,': VK MASK ',24i1)
737     if(ima.ne.0.and.verbose.eq.1)write(*,100)eventn(1),iv
738     $ ,(mask_vk_ev(iv,ik),ik=1,nva1_view)
739    
740     return
741     end
742    
743     *---***---***---***---***---***---***---***---***
744     *
745     *
746     *
747     *
748     *
749     *---***---***---***---***---***---***---***---***
750     subroutine subtractped(iv)
751    
752     include 'commontracker.f'
753     include 'level1.f'
754     include 'calib.f'
755     include 'common_reduction.f'
756 mocchiut 1.1
757 pam-fi 1.25 do is=1,nstrips_view !loop on strips (1)
758     if(mod(iv,2).eq.1) then
759     C=== > Y view
760     c print*,iv,nvk(is),nst(is),adc(iv,nvk(is),nst(is))
761     c $ ,cn(iv,nvk(is))
762     c $ ,pedestal(iv,nvk(is),nst(is))
763     value(is)= -(DBLE(adc(iv,nvk(is),nst(is)))
764     $ -cn(iv,nvk(is))-pedestal(iv,nvk(is),nst(is)))
765     $ *mask(iv,nvk(is),nst(is))
766     clseedcut(is)=clcuty*sigma(iv,nvk(is),nst(is))
767     $ *mask(iv,nvk(is),nst(is))
768     clinclcut(is)=incuty*sigma(iv,nvk(is),nst(is))
769     $ *mask(iv,nvk(is),nst(is))
770     sat(is)=0
771 pam-fi 1.26 if( adc(iv,nvk(is),nst(is)).lt.adc_saty )
772     $ sat(is)=mask(iv,nvk(is),nst(is))
773 pam-fi 1.25 else
774     C=== > X view
775     value(is)= (DBLE(adc(iv,nvk(is),nst(is)))
776     $ -cn(iv,nvk(is))-pedestal(iv,nvk(is),nst(is)))
777     $ *mask(iv,nvk(is),nst(is))
778     clseedcut(is)=clcutx*sigma(iv,nvk(is),nst(is))
779     $ *mask(iv,nvk(is),nst(is))
780     clinclcut(is)=incutx*sigma(iv,nvk(is),nst(is))
781     $ *mask(iv,nvk(is),nst(is))
782     sat(is)=0
783 pam-fi 1.26 if( adc(iv,nvk(is),nst(is)).gt.adc_satx )
784     $ sat(is)=mask(iv,nvk(is),nst(is))
785 pam-fi 1.25 endif
786     enddo !end loop on strips (1)
787    
788    
789     return
790     end
791     *---***---***---***---***---***---***---***---***
792     *
793     *
794     *
795     *
796     *
797     *---***---***---***---***---***---***---***---***
798 pam-fi 1.13 c$$$ subroutine stripmask
799     c$$$
800     c$$$* this routine set va1 and single-strip masks,
801     c$$$* on the basis of the VA1 mask saved in the DB
802     c$$$*
803     c$$$* mask(nviews,nva1_view,nstrips_va1) !strip mask
804     c$$$* mask_vk(nviews,nva1_view) !VA1 mask
805     c$$$*
806     c$$$ include 'commontracker.f'
807     c$$$ include 'level1.f'
808     c$$$ include 'common_reduction.f'
809     c$$$ include 'calib.f'
810     c$$$
811     c$$$* init mask
812     c$$$ do iv=1,nviews
813     c$$$ do ivk=1,nva1_view
814     c$$$ do is=1,nstrips_va1
815     c$$$c mask(iv,ivk,is) = mask_vk(iv,ivk)
816     c$$$ if( mask_vk(iv,ivk) .ne. -1)then
817     c$$$ mask(iv,ivk,is) = 1
818     c$$$ $ * mask_vk(iv,ivk) !from DB
819     c$$$ $ * mask_vk_ev(iv,ivk) !from <SIG>
820     c$$$ $ * mask_vk_run(iv,ivk) !from CN
821     c$$$ else
822     c$$$ mask(iv,ivk,is) = -1
823     c$$$ $ * mask_vk(iv,ivk) !from DB
824     c$$$ $ * mask_vk_ev(iv,ivk) !from CN
825     c$$$ endif
826     c$$$ enddo
827     c$$$ enddo
828     c$$$ enddo
829     c$$$
830     c$$$
831     c$$$ return
832     c$$$ end
833    
834     subroutine stripmask(iv,ivk)
835 mocchiut 1.1
836 pam-fi 1.18 * -----------------------------------------------
837 mocchiut 1.1 * this routine set va1 and single-strip masks,
838     * on the basis of the VA1 mask saved in the DB
839     *
840     * mask(nviews,nva1_view,nstrips_va1) !strip mask
841     * mask_vk(nviews,nva1_view) !VA1 mask
842 pam-fi 1.18 * -----------------------------------------------
843 mocchiut 1.1 include 'commontracker.f'
844 pam-fi 1.5 include 'level1.f'
845 pam-fi 1.4 include 'common_reduction.f'
846 mocchiut 1.1 include 'calib.f'
847    
848     * init mask
849 pam-fi 1.13 do is=1,nstrips_va1
850 pam-fi 1.18 * --------------------------------------------------------
851     * if VA1-mask from DB is 0 or 1, three masks are combined:
852     * - from DB (a-priori mask)
853     * - run-based (chip declared bad on the basis of <SIG>)
854     * - event-based (failure in CN computation)
855     * --------------------------------------------------------
856 pam-fi 1.20 c print*,iv,ivk
857     c $ ,mask_vk(iv,ivk),mask_vk_ev(iv,ivk),mask_vk_run(iv,ivk)
858 pam-fi 1.13 if( mask_vk(iv,ivk) .ne. -1)then
859     mask(iv,ivk,is) = 1
860 pam-fi 1.18 $ * mask_vk(iv,ivk) !from DB
861     $ * mask_vk_ev(iv,ivk) !from <SIG>
862 pam-fi 1.13 $ * mask_vk_run(iv,ivk) !from CN
863 pam-fi 1.18 * -----------------------------------------------------------
864     * if VA1-mask from DB is -1 only event-based mask is applied
865     * -----------------------------------------------------------
866 pam-fi 1.13 else
867     mask(iv,ivk,is) = -1
868 pam-fi 1.18 $ * mask_vk(iv,ivk) !from DB
869     $ * mask_vk_ev(iv,ivk) !from CN
870 pam-fi 1.13 endif
871 mocchiut 1.1 enddo
872 pam-fi 1.13
873    
874 mocchiut 1.1 return
875     end

  ViewVC Help
Powered by ViewVC 1.1.23