/[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.25 - (hide annotations) (download)
Wed Oct 22 15:17:40 2008 UTC (16 years, 1 month ago) by pam-fi
Branch: MAIN
Changes since 1.24: +167 -77 lines
fixed bug in TrkHough track + new method UnpackError()

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     if(debug.eq.1)print*,'view ',iv,' #clusters ', nclstr_view
247     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.10 itemp = is
397     fsat = 0 ! first saturated strip
398     lsat = 0 ! last saturated strip
399 mocchiut 1.1 if(itemp.eq.last) goto 230 !estremo...
400 pam-fi 1.10 c ------------------------
401     c search for first maximum
402     c ------------------------
403 pam-fi 1.6 do while(
404     $ value(itemp).le.value(itemp+1)
405     $ .and.value(itemp+1).gt.clseedcut(itemp+1))
406 pam-fi 1.10 itemp = itemp+1
407     if(itemp.eq.last) goto 230 !stops if reaches last strip
408     if(sat(itemp).eq.1) goto 230 !stop if reaches a saturated strip
409 mocchiut 1.1 enddo ! of the ladder
410     230 continue
411 pam-fi 1.10 c -----------------------------
412     c check if strips are saturated
413     c -----------------------------
414     if( sat(itemp).eq.1 )then
415     fsat = itemp
416     lsat = itemp
417     if(itemp.eq.last) goto 231 !estremo...
418     do while( sat(itemp+1).eq.1 )
419     itemp = itemp+1
420     lsat = itemp
421     if(itemp.eq.last) goto 231 !stops if reaches last strip
422     enddo
423     endif
424     231 continue
425     c---------------------------------------------------------------------------
426 mocchiut 1.1 c fownd SEED!!!
427 pam-fi 1.10 c (if there are saturated strips, the cluster is centered in the middle)
428     c---------------------------------------------------------------------------
429     if(fsat.eq.0.and.lsat.eq.0)then
430     iseed = itemp ! <<< SEED
431     else
432     iseed = int((lsat+fsat)/2) ! <<< SEED
433     c$$$ print*,'saturated strips ',fsat,lsat,iseed
434     c$$$ print*,'--> ',(value(ii),ii=fsat,lsat)
435     endif
436     c---------------------------------------------------------------
437 mocchiut 1.1 c after finding a cluster seed, checks also adjacent strips,
438 pam-fi 1.10 C and tags the ones exceeding clinclcut
439     c---------------------------------------------------------------
440 mocchiut 1.1 ir=iseed !indici destro
441     il=iseed ! e sinistro
442    
443     rmax=ir !estremo destro del cluster
444     lmax=il ! e sinistro
445    
446     rstop=0 !initialize flags used to exit from
447     lstop=0 ! inclusion loop
448    
449     do while(lstop.eq.0.or.rstop.eq.0) !shifts left and right from
450 pam-fi 1.10
451    
452     ir=ir+1 !index for right side
453     il=il-1 !index for left side
454 mocchiut 1.1 c------------------------------------------------------------------------
455     c checks for last or first strip of the ladder
456     c------------------------------------------------------------------------
457 pam-fi 1.10 if( ir.gt.last ) rstop = 1
458     if( il.lt.first ) lstop = 1
459 mocchiut 1.1
460     c------------------------------------------------------------------------
461 pam-fi 1.10 c add strips exceeding inclusion cut
462 mocchiut 1.1 c------------------------------------------------------------------------
463 pam-fi 1.10 if( (rmax-lmax+1).ge.nclstrp )goto 210 !exits inclusion loop
464    
465 pam-fi 1.21 if(rstop.eq.0) then !if right cluster border has not been reached
466 pam-fi 1.10 if(value(ir).gt.clinclcut(ir)) then
467     rmax=ir !include a strip on the right
468 mocchiut 1.1 else
469 pam-fi 1.10 rstop=1 !cluster right end
470     endif
471 mocchiut 1.1 endif
472 pam-fi 1.10
473     if( (rmax-lmax+1).ge.nclstrp )goto 210 !exits inclusion loop
474    
475 pam-fi 1.21 if(lstop.eq.0) then !if left cluster border has not been reached
476 mocchiut 1.1 if(value(il).gt.clinclcut(il)) then
477 pam-fi 1.10 lmax=il !include a strip on the left
478 mocchiut 1.1 else
479 pam-fi 1.10 lstop=1 !cluster left end
480 mocchiut 1.1 endif
481     endif
482    
483 pam-fi 1.21 c if( (rmax-lmax+1).ge.nclstrp )goto 210 !exits inclusion loop
484    
485 mocchiut 1.1 enddo !ends strip inclusion loop
486 pam-fi 1.10 goto 211
487 mocchiut 1.1 210 continue !jumps here if more than nclstrp have been included
488 pam-fi 1.10 c print*,'>>> nclstrp! '
489     211 continue
490     c-----------------------------------------
491     c end of inclusion loop!
492     c-----------------------------------------
493    
494     c------------------------------------------------------------------------
495     c adjust the cluster in order to have at least a strip around the seed
496     c------------------------------------------------------------------------
497     if(iseed.eq.lmax.and.lmax.ne.first)then
498     lmax = lmax-1
499     if( (rmax-lmax+1).gt.nclstrp )rmax=rmax-1
500     endif
501     if(iseed.eq.rmax.and.rmax.ne.last )then
502     rmax = rmax+1
503     if( (rmax-lmax+1).gt.nclstrp )lmax=lmax+1
504 mocchiut 1.1 endif
505 pam-fi 1.21 c-------------------------------------------------------------------------------
506     c adjust the cluster in order to have at least ANOTHER strip around the seed
507     c-------------------------------------------------------------------------------
508 pam-fi 1.24 if(iseed-1.eq.lmax.and.lmax.ne.first)then
509     lmax = lmax-1
510     if( (rmax-lmax+1).gt.nclstrp )rmax=rmax-1
511     endif
512     if(iseed+1.eq.rmax.and.rmax.ne.last )then
513     rmax = rmax+1
514     if( (rmax-lmax+1).gt.nclstrp )lmax=lmax+1
515     endif
516 pam-fi 1.21 c---------------------------------------------------
517     c now we have 5 stored-strips around the maximum
518     c---------------------------------------------------
519 pam-fi 1.10
520 mocchiut 1.1 c------------------------------------------------------------------------
521 pam-fi 1.10 c adjust the cluster in order to store a minimum number of strips
522     c------------------------------------------------------------------------
523     do while( (rmax-lmax+1).lt.nclstrpmin )
524    
525     vl = -99999
526     vr = -99999
527     if(lmax-1.ge.first) vl = value(lmax-1)
528     if(rmax+1.le.last ) vr = value(rmax+1)
529     if(vr.ge.vl) then
530     rmax = rmax+1
531     else
532     lmax = lmax-1
533 mocchiut 1.1 endif
534 pam-fi 1.10
535     enddo
536 mocchiut 1.1
537     c--------------------------------------------------------
538 pam-fi 1.10 c store cluster info
539 mocchiut 1.1 c--------------------------------------------------------
540 pam-fi 1.5 nclstr_view = nclstr_view + 1 !cluster number
541 pam-fi 1.10
542 pam-fi 1.5 if(nclstr_view.gt.nclstrmax_view) then !too many clusters for the view:
543 pam-fi 1.10 c$$$ if(verbose) print*,'Event ',eventn(1),
544     c$$$ $ ': more than ',nclstrmax_view
545     c$$$ $ ,' clusters on view ',iv
546 mocchiut 1.1 flag_shower = .true.
547     goto 2000
548     endif
549 pam-fi 1.5
550 pam-fi 1.10 ladder_view(nclstr_view) = nld(iseed,iv)
551     maxs_view(nclstr_view) = iseed
552 pam-fi 1.5 rmax_view(nclstr_view) = rmax
553     lmax_view(nclstr_view) = lmax
554 pam-fi 1.23 c mult_view(nclstr_view) = rmax-lmax+1
555     mult_view(nclstr_view) = 0
556     do ii=lmax,rmax
557     if(value(ii).gt.clinclcut(ii))
558     $ mult_view(nclstr_view) = mult_view(nclstr_view)+1
559     enddo
560    
561 pam-fi 1.5
562 pam-fi 1.10 c$$$ if(rmax-lmax+1.gt.25)
563     c$$$ $ print*,'view ',iv
564     c$$$ $ ,' cl ',nclstr_view,' mult ',rmax-lmax+1
565     c------------------------------------------------------------------------
566 pam-fi 1.11 c search for a double peak inside the cluster
567 pam-fi 1.10 c------------------------------------------------------------------------
568     inext = rmax+1 !<< index where to start new-cluster search
569    
570     vmax = 0
571     vmin = value(iseed)
572     imax = iseed
573     imin = iseed
574     do iss = max(iseed+1,lsat+1),rmax
575     if( value(iss).lt.vmin )then
576     if( imax.ne.iseed )goto 221 !found dowble peek
577     imin = iss
578     vmin = value(iss)
579     else
580     delta = value(iss) - value(imin)
581     cut = sqrt(clinclcut(iss)**2 + clinclcut(imin)**2)
582     if(
583     $ delta.gt.cut .and.
584     $ value(iss).gt.clseedcut(iss).and.
585     $ .true.)then
586     if( value(iss).gt.vmax )then
587     imax = iss
588     vmax = value(iss)
589     else
590     goto 221 !found dowble peek
591     endif
592     endif
593     endif
594     enddo
595     221 continue
596    
597     if(imax.gt.iseed)then
598     inext = imax !<< index where to start new-cluster search
599     c$$$ print*,'--- double peek ---'
600     c$$$ print*,(value(ii),ii=lmax,rmax)
601     c$$$ print*,'seed ',iseed,' imin ',imin,' imax ',imax
602     endif
603 mocchiut 1.1 c--------------------------------------------------------
604 pam-fi 1.2 c
605 mocchiut 1.1 c--------------------------------------------------------
606     endif !end possible seed conditio
607     220 continue !jumps here to skip strips left of last seed
608    
609     enddo ! end loop on strips
610     enddo !end loop on ladders
611     2000 continue
612     return
613     end
614    
615    
616     *---***---***---***---***---***---***---***---***
617     *
618     *
619     *
620     *
621     *
622     *---***---***---***---***---***---***---***---***
623    
624 pam-fi 1.25 subroutine savecluster(iv)
625 pam-fi 1.5 *
626     * (080/2006 Elena Vannuccini)
627     * Save the clusters view by view
628    
629     include 'commontracker.f'
630     include 'level1.f'
631     include 'calib.f'
632     include 'common_reduction.f'
633    
634     integer CLlength !lunghezza in strip del cluster
635    
636     do ic=1,nclstr_view
637    
638     nclstr1 = nclstr1+1
639     view(nclstr1) = iv
640     ladder(nclstr1) = ladder_view(ic)
641     maxs(nclstr1) = maxs_view(ic)
642     mult(nclstr1) = mult_view(ic)
643    
644     c posizione dell'inizio del cluster nell' array clsignal
645     indstart(nclstr1) = ind
646     c posizione del cluster seed nell'array clsignal
647     indmax(nclstr1) = indstart(nclstr1)
648     $ +( maxs_view(ic) - lmax_view(ic) )
649    
650     CLlength = rmax_view(ic) - lmax_view(ic) + 1 !numero di strip salvate
651     totCLlength = totCLlength + CLlength
652 pam-fi 1.16 sgnl(nclstr1) = 0
653 pam-fi 1.5 do j=lmax_view(ic),rmax_view(ic) !stores sequentially cluter strip values in
654    
655     clsignal(ind) = value(j) ! clsignal array
656 pam-fi 1.20 c$$$ print*,ind,clsignal(ind)
657 pam-fi 1.5 ivk=nvk(j)
658     ist=nst(j)
659    
660     clsigma(ind) = sigma(iv,ivk,ist)
661     cladc(ind) = adc(iv,ivk,ist)
662     clbad(ind) = bad(iv,ivk,ist)
663     c clped(ind) = pedestal(iv,ivk,ist)
664    
665     ind=ind+1
666     c if(value(j).gt.0)
667     if(value(j).gt.clinclcut(j))
668 pam-fi 1.16 $ sgnl(nclstr1) = sgnl(nclstr1) + value(j) !cluster charge
669 pam-fi 1.5 enddo
670    
671 pam-fi 1.25 c$$$ print*,'view ',iv,' -- savecluster -- nclstr1: '
672 pam-fi 1.20 c$$$ $ ,nclstr1,maxs(nclstr1),mult(nclstr1),sgnl(nclstr1)
673     c$$$ print*,'----------------------'
674    
675 pam-fi 1.5 enddo
676    
677     return
678     end
679     *---***---***---***---***---***---***---***---***
680     *
681     *
682     *
683     *
684     *
685     *---***---***---***---***---***---***---***---***
686    
687 pam-fi 1.25 subroutine evaluatecn(iv)
688    
689     include 'commontracker.f'
690     include 'level0.f'
691     include 'level1.f'
692     include 'common_reduction.f'
693     include 'calib.f'
694    
695     ima=0
696     do ik=1,nva1_view
697     cn(iv,ik) = 0
698     cnrms(iv,ik) = 0
699     cnn(iv,ik) = -1
700     iflag = 0
701     mask_vk_ev(iv,ik) = 1
702     call stripmask(iv,ik) !compute mask(i,j,k), combining VA1-masks
703     * --------------------------------------
704     * if chip is not masked ---> evaluate CN
705     * --------------------------------------
706     if( mask(iv,ik,1).eq.1 ) then !!!NBNB mask per la striscia 1 !!!!!!!!
707     call cncomp(iv,ik,iflag)
708     if(iflag.ne.0)then
709     ima=ima+1
710     mask_vk_ev(iv,ik)=0
711     ierror = 220
712     endif
713     call stripmask(iv,ik) !compute mask(i,j,k), combining VA1-masks
714     endif
715     enddo
716     100 format(' * WARNING * Event ',i7,' view',i3,': VK MASK ',24i1)
717     if(ima.ne.0.and.verbose.eq.1)write(*,100)eventn(1),iv
718     $ ,(mask_vk_ev(iv,ik),ik=1,nva1_view)
719    
720     return
721     end
722    
723     *---***---***---***---***---***---***---***---***
724     *
725     *
726     *
727     *
728     *
729     *---***---***---***---***---***---***---***---***
730     subroutine subtractped(iv)
731    
732     include 'commontracker.f'
733     include 'level1.f'
734     include 'calib.f'
735     include 'common_reduction.f'
736 mocchiut 1.1
737 pam-fi 1.25 do is=1,nstrips_view !loop on strips (1)
738     if(mod(iv,2).eq.1) then
739     C=== > Y view
740     c print*,iv,nvk(is),nst(is),adc(iv,nvk(is),nst(is))
741     c $ ,cn(iv,nvk(is))
742     c $ ,pedestal(iv,nvk(is),nst(is))
743     value(is)= -(DBLE(adc(iv,nvk(is),nst(is)))
744     $ -cn(iv,nvk(is))-pedestal(iv,nvk(is),nst(is)))
745     $ *mask(iv,nvk(is),nst(is))
746     clseedcut(is)=clcuty*sigma(iv,nvk(is),nst(is))
747     $ *mask(iv,nvk(is),nst(is))
748     clinclcut(is)=incuty*sigma(iv,nvk(is),nst(is))
749     $ *mask(iv,nvk(is),nst(is))
750     sat(is)=0
751     if( adc(iv,nvk(is),nst(is)).lt.adc_saty )sat(is)=1
752     else
753     C=== > X view
754     value(is)= (DBLE(adc(iv,nvk(is),nst(is)))
755     $ -cn(iv,nvk(is))-pedestal(iv,nvk(is),nst(is)))
756     $ *mask(iv,nvk(is),nst(is))
757     clseedcut(is)=clcutx*sigma(iv,nvk(is),nst(is))
758     $ *mask(iv,nvk(is),nst(is))
759     clinclcut(is)=incutx*sigma(iv,nvk(is),nst(is))
760     $ *mask(iv,nvk(is),nst(is))
761     sat(is)=0
762     if( adc(iv,nvk(is),nst(is)).gt.adc_satx )sat(is)=1
763     endif
764     enddo !end loop on strips (1)
765    
766    
767     return
768     end
769     *---***---***---***---***---***---***---***---***
770     *
771     *
772     *
773     *
774     *
775     *---***---***---***---***---***---***---***---***
776 pam-fi 1.13 c$$$ subroutine stripmask
777     c$$$
778     c$$$* this routine set va1 and single-strip masks,
779     c$$$* on the basis of the VA1 mask saved in the DB
780     c$$$*
781     c$$$* mask(nviews,nva1_view,nstrips_va1) !strip mask
782     c$$$* mask_vk(nviews,nva1_view) !VA1 mask
783     c$$$*
784     c$$$ include 'commontracker.f'
785     c$$$ include 'level1.f'
786     c$$$ include 'common_reduction.f'
787     c$$$ include 'calib.f'
788     c$$$
789     c$$$* init mask
790     c$$$ do iv=1,nviews
791     c$$$ do ivk=1,nva1_view
792     c$$$ do is=1,nstrips_va1
793     c$$$c mask(iv,ivk,is) = mask_vk(iv,ivk)
794     c$$$ if( mask_vk(iv,ivk) .ne. -1)then
795     c$$$ mask(iv,ivk,is) = 1
796     c$$$ $ * mask_vk(iv,ivk) !from DB
797     c$$$ $ * mask_vk_ev(iv,ivk) !from <SIG>
798     c$$$ $ * mask_vk_run(iv,ivk) !from CN
799     c$$$ else
800     c$$$ mask(iv,ivk,is) = -1
801     c$$$ $ * mask_vk(iv,ivk) !from DB
802     c$$$ $ * mask_vk_ev(iv,ivk) !from CN
803     c$$$ endif
804     c$$$ enddo
805     c$$$ enddo
806     c$$$ enddo
807     c$$$
808     c$$$
809     c$$$ return
810     c$$$ end
811    
812     subroutine stripmask(iv,ivk)
813 mocchiut 1.1
814 pam-fi 1.18 * -----------------------------------------------
815 mocchiut 1.1 * this routine set va1 and single-strip masks,
816     * on the basis of the VA1 mask saved in the DB
817     *
818     * mask(nviews,nva1_view,nstrips_va1) !strip mask
819     * mask_vk(nviews,nva1_view) !VA1 mask
820 pam-fi 1.18 * -----------------------------------------------
821 mocchiut 1.1 include 'commontracker.f'
822 pam-fi 1.5 include 'level1.f'
823 pam-fi 1.4 include 'common_reduction.f'
824 mocchiut 1.1 include 'calib.f'
825    
826     * init mask
827 pam-fi 1.13 do is=1,nstrips_va1
828 pam-fi 1.18 * --------------------------------------------------------
829     * if VA1-mask from DB is 0 or 1, three masks are combined:
830     * - from DB (a-priori mask)
831     * - run-based (chip declared bad on the basis of <SIG>)
832     * - event-based (failure in CN computation)
833     * --------------------------------------------------------
834 pam-fi 1.20 c print*,iv,ivk
835     c $ ,mask_vk(iv,ivk),mask_vk_ev(iv,ivk),mask_vk_run(iv,ivk)
836 pam-fi 1.13 if( mask_vk(iv,ivk) .ne. -1)then
837     mask(iv,ivk,is) = 1
838 pam-fi 1.18 $ * mask_vk(iv,ivk) !from DB
839     $ * mask_vk_ev(iv,ivk) !from <SIG>
840 pam-fi 1.13 $ * mask_vk_run(iv,ivk) !from CN
841 pam-fi 1.18 * -----------------------------------------------------------
842     * if VA1-mask from DB is -1 only event-based mask is applied
843     * -----------------------------------------------------------
844 pam-fi 1.13 else
845     mask(iv,ivk,is) = -1
846 pam-fi 1.18 $ * mask_vk(iv,ivk) !from DB
847     $ * mask_vk_ev(iv,ivk) !from CN
848 pam-fi 1.13 endif
849 mocchiut 1.1 enddo
850 pam-fi 1.13
851    
852 mocchiut 1.1 return
853     end

  ViewVC Help
Powered by ViewVC 1.1.23