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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.18 - (show annotations) (download)
Fri Apr 27 10:39:58 2007 UTC (17 years, 7 months ago) by pam-fi
Branch: MAIN
CVS Tags: v3r04, v3r03
Changes since 1.17: +33 -18 lines
v3r00: new hough parameters, new variables, and other things...

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 subroutine reductionflight(ierror)
14
15 include 'commontracker.f'
16 include 'level0.f'
17 include 'level1.f'
18 include 'common_reduction.f'
19 include 'calib.f'
20
21 data eventn_old/nviews*0/
22
23 integer ierror
24 ierror = 0
25
26 c$$$ debug = .true.
27 c$$$ verbose = .true.
28 c$$$ warning = .true.
29
30 call init_level1
31
32 if(debug)print*,'-- check LEVEL0 status'
33
34 c good1 = good0
35 c--------------------------------------------------
36 c check the LEVEL0 event status for missing
37 c sections or DSP alarms
38 c ==> set the variable GOOD1(12)
39 c--------------------------------------------------
40 do iv=1,nviews
41 if(DSPnumber(iv).gt.0.and.DSPnumber(iv).le.12)then
42 c ------------------------
43 c GOOD
44 c ------------------------
45 GOOD1(DSPnumber(iv))=0 !OK
46 c ------------------------
47 c CRC error
48 c ------------------------
49 if(crc(iv).eq.1) then
50 GOOD1(DSPnumber(iv)) = 2
51 goto 18 !next view
52 endif
53 c ------------------------
54 c online-software alarm
55 c ------------------------
56 if(
57 $ fl1(iv).ne.0.or.
58 $ fl2(iv).ne.0.or.
59 $ fl3(iv).ne.0.or.
60 $ fl4(iv).ne.0.or.
61 $ fl5(iv).ne.0.or.
62 $ fl6(iv).ne.0.or.
63 $ fc(iv).ne.0.or.
64 $ DATAlength(iv).eq.0.or.
65 $ .false.)then
66 GOOD1(DSPnumber(iv))=3
67 goto 18
68 endif
69 c ------------------------
70 c DSP-counter jump
71 c ------------------------
72 if(
73 $ eventn_old(iv).ne.0.and. !first event in this file
74 $ eventn(iv).ne.1.and. !first event in run
75 $ good_old(DSPnumber(iv)).ne.0.and. !previous event corrupted
76 $ .true.)then
77
78 if(eventn(iv).ne.(eventn_old(iv)+1))then
79 GOOD1(DSPnumber(iv))=4
80 goto 18
81 endif
82
83 endif
84 c ------------------------
85 18 continue
86 endif
87 enddo
88
89 ngood = 0
90 do iv = 1,nviews
91 eventn_old(iv) = eventn(iv)
92 good_old(iv) = good1(iv)
93 ngood = ngood + good1(iv)
94 enddo
95 if(debug.and.ngood.ne.0)print*,'* WARNING * LEVEL0 event status: '
96 $ ,(good1(i),i=1,nviews)
97 c--------------------------------------------------
98 c read the variable DATATRACKER from LEVEL0
99 c and fill the variable ADC (invertin view 11)
100 c--------------------------------------------------
101
102 if(debug)print*,'-- fill ADC vectors'
103
104 call filladc(iflag)
105 if(iflag.ne.0)then
106 ierror = 220
107 endif
108
109 c--------------------------------------------------
110 c computes common noise for each VA1
111 c (excluding strips with signal,
112 c tagged with the flag CLSTR)
113 c--------------------------------------------------
114 if(debug)print*,'-- compute CN'
115
116 do iv=1,nviews
117 ima=0
118 do ik=1,nva1_view
119 cn(iv,ik) = 0
120 cnrms(iv,ik) = 0
121 cnn(iv,ik) = -1
122 iflag = 0
123 mask_vk_ev(iv,ik) = 1
124 call stripmask(iv,ik) !compute mask(i,j,k), combining VA1-masks
125 * --------------------------------------
126 * if chip is not masked ---> evaluate CN
127 * --------------------------------------
128 if( mask(iv,ik,1).eq.1 ) then !!!NBNB mask per la striscia 1 !!!!!!!!
129 call cncomp(iv,ik,iflag)
130 if(iflag.ne.0)then
131 ima=ima+1
132 mask_vk_ev(iv,ik)=0
133 ierror = 220
134 endif
135 call stripmask(iv,ik) !compute mask(i,j,k), combining VA1-masks
136 endif
137 enddo
138 100 format(' * WARNING * Event ',i7,' view',i3,': VK MASK ',24i1)
139 if(ima.ne.0.and.debug)write(*,100)eventn(1),iv
140 $ ,(mask_vk_ev(iv,ik),ik=1,nva1_view)
141 c if(ima.ne.0)write(*,100)eventn(1),iv
142 c $ ,(mask_vk_ev(iv,ik),ik=1,nva1_view)
143 enddo
144
145 cc call stripmask !compute mask(i,j,k), combining mask_vk_ev and mask_vk
146
147 c---------------------------------------------
148 c loops on views, VA1 and strips,
149 c and computes strips signals using
150 c badstrip, pedestals, and
151 c sigma informations from histograms
152 c---------------------------------------------
153 ind=1 !clsignal array index
154
155 if(debug)print*,'-- search clusters'
156 do iv=1,nviews !loop on views
157 do is=1,nstrips_view !loop on strips (1)
158 if(mod(iv,2).eq.1) then
159 C=== > Y view
160 value(is)= -(DBLE(adc(iv,nvk(is),nst(is)))
161 $ -cn(iv,nvk(is))-pedestal(iv,nvk(is),nst(is)))
162 $ *mask(iv,nvk(is),nst(is))
163 clseedcut(is)=clcuty*sigma(iv,nvk(is),nst(is))
164 $ *mask(iv,nvk(is),nst(is))
165 clinclcut(is)=incuty*sigma(iv,nvk(is),nst(is))
166 $ *mask(iv,nvk(is),nst(is))
167 sat(is)=0
168 if( adc(iv,nvk(is),nst(is)).lt.adc_saty )sat(is)=1
169 else
170 C=== > X view
171 value(is)= (DBLE(adc(iv,nvk(is),nst(is)))
172 $ -cn(iv,nvk(is))-pedestal(iv,nvk(is),nst(is)))
173 $ *mask(iv,nvk(is),nst(is))
174 clseedcut(is)=clcutx*sigma(iv,nvk(is),nst(is))
175 $ *mask(iv,nvk(is),nst(is))
176 clinclcut(is)=incutx*sigma(iv,nvk(is),nst(is))
177 $ *mask(iv,nvk(is),nst(is))
178 sat(is)=0
179 if( adc(iv,nvk(is),nst(is)).gt.adc_satx )sat(is)=1
180 endif
181 enddo !end loop on strips (1)
182 call search_cluster(iv)
183
184 if(.not.flag_shower)then
185 call save_cluster(iv)
186 if(debug)print*,'view ',iv,' #clusters ', nclstr_view
187 else
188 fshower(iv) = 1
189 c GOOD1(DSPnumber(iv)) = 11 !AHAHAHAHA ORRORE!!
190 GOOD1(iv) = 11
191 101 format(' * WARNING * Event ',i7,' view',i3
192 $ ,' #clusters > ',i5,' --> MASKED')
193 if(debug)write(*,101)eventn(1),iv,nclstrmax_view
194 endif
195 enddo ! end loop on views
196 do iv=1,nviews
197 do ik=1,nva1_view
198 cnev(iv,ik) = cn(iv,ik) !assigns computed CN to ntuple variables
199 cnrmsev(iv,ik) = cnrms(iv,ik) !assigns computed CN to ntuple variables
200 cnnev(iv,ik) = cnn(iv,ik) !assigns computed CN to ntuple variables
201 enddo
202 enddo
203 C---------------------------------------------
204 C come here if GOOD1=0
205 C or the event has too many clusters
206 C---------------------------------------------
207 200 continue
208
209 ngood = 0
210 do iv = 1,nviews
211 ngood = ngood + good1(iv)
212 enddo
213 if(debug.and.ngood.ne.0)print*,'* WARNING * Event ',eventn(1)
214 $ ,':LEVEL1 event status: '
215 $ ,(good1(i),i=1,nviews)
216 c------------------------------------------------------------------------
217 c
218 c closes files and exits
219 c
220 c------------------------------------------------------------------------
221 RETURN
222 END
223
224 ***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***
225 *
226 *
227 *
228 *
229 *
230 *
231 *
232 *
233 *
234 ***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***
235
236
237 subroutine init_level1
238
239 include 'commontracker.f'
240 include 'level1.f'
241 include 'level0.f'
242
243 c good1 = 0
244 do iv=1,12
245 good1(iv) = 1 !missing packet
246 enddo
247 nclstr1 = 0
248 totCLlength = 0
249 do ic=1,nclstrmax
250 view(ic) = 0
251 ladder(ic) = 0
252 indstart(ic) = 0
253 indmax(ic) = 0
254 maxs(ic) = 0
255 mult(ic) = 0
256 sgnl(ic) = 0
257 whichtrack(ic) = 0 !assigned @ level2
258
259 enddo
260 do id=1,maxlength !???
261 clsignal(id) = 0.
262 clsigma(id) = 0.
263 cladc(id) = 0.
264 clbad(id) = 0.
265 enddo
266 do iv=1,nviews
267 c crc1(iv)=0
268 do ik=1,nva1_view
269 cnev(iv,ik) = 0
270 cnnev(iv,ik) = 0
271 enddo
272 fshower(iv) = 0
273 enddo
274
275 return
276 end
277
278 *---***---***---***---***---***---***---***---***
279 *
280 *
281 *
282 *
283 *
284 *---***---***---***---***---***---***---***---***
285
286 subroutine search_cluster(iv)
287
288 include 'commontracker.f'
289 include 'level0.f'
290 include 'level1.f'
291 include 'calib.f'
292
293 include 'common_reduction.f'
294
295
296 c local variables
297 integer rmax,lmax !estremi del cluster
298 integer rstop,lstop
299 integer first,last
300 integer fsat,lsat
301
302 external nst
303
304 iseed=-999 !cluster seed index initialization
305
306 inext=-999 !index where to start new cluster search
307
308 flag_shower = .false.
309 nclstr_view=0
310
311 do jl=1,nladders_view !1..3 !loops on ladders
312
313 first = 1+nstrips_ladder*(jl-1) !1,1025,2049
314 last = nstrips_ladder*jl !1024,2048,3072
315
316 * X views have 1018 strips instead of 1024
317 if(mod(iv,2).eq.0) then
318 first=first+3
319 last=last-3
320 endif
321
322 do is=first,last !loop on strips in each ladder
323
324 c---------------------------------------------
325 c new-cluster search starts at index inext
326 c---------------------------------------------
327 if(is.lt.inext) goto 220 ! next strip
328
329 if(value(is).gt.clseedcut(is)) then
330 c-----------------------------------------
331 c possible SEED...
332 c-----------------------------------------
333 itemp = is
334 fsat = 0 ! first saturated strip
335 lsat = 0 ! last saturated strip
336 if(itemp.eq.last) goto 230 !estremo...
337 c ------------------------
338 c search for first maximum
339 c ------------------------
340 do while(
341 $ value(itemp).le.value(itemp+1)
342 $ .and.value(itemp+1).gt.clseedcut(itemp+1))
343 itemp = itemp+1
344 if(itemp.eq.last) goto 230 !stops if reaches last strip
345 if(sat(itemp).eq.1) goto 230 !stop if reaches a saturated strip
346 enddo ! of the ladder
347 230 continue
348 c -----------------------------
349 c check if strips are saturated
350 c -----------------------------
351 if( sat(itemp).eq.1 )then
352 fsat = itemp
353 lsat = itemp
354 if(itemp.eq.last) goto 231 !estremo...
355 do while( sat(itemp+1).eq.1 )
356 itemp = itemp+1
357 lsat = itemp
358 if(itemp.eq.last) goto 231 !stops if reaches last strip
359 enddo
360 endif
361 231 continue
362 c---------------------------------------------------------------------------
363 c fownd SEED!!!
364 c (if there are saturated strips, the cluster is centered in the middle)
365 c---------------------------------------------------------------------------
366 if(fsat.eq.0.and.lsat.eq.0)then
367 iseed = itemp ! <<< SEED
368 else
369 iseed = int((lsat+fsat)/2) ! <<< SEED
370 c$$$ print*,'saturated strips ',fsat,lsat,iseed
371 c$$$ print*,'--> ',(value(ii),ii=fsat,lsat)
372 endif
373 c---------------------------------------------------------------
374 c after finding a cluster seed, checks also adjacent strips,
375 C and tags the ones exceeding clinclcut
376 c---------------------------------------------------------------
377 ir=iseed !indici destro
378 il=iseed ! e sinistro
379
380 rmax=ir !estremo destro del cluster
381 lmax=il ! e sinistro
382
383 rstop=0 !initialize flags used to exit from
384 lstop=0 ! inclusion loop
385
386 do while(lstop.eq.0.or.rstop.eq.0) !shifts left and right from
387
388
389 ir=ir+1 !index for right side
390 il=il-1 !index for left side
391 c------------------------------------------------------------------------
392 c checks for last or first strip of the ladder
393 c------------------------------------------------------------------------
394 if( ir.gt.last ) rstop = 1
395 if( il.lt.first ) lstop = 1
396
397 c------------------------------------------------------------------------
398 c add strips exceeding inclusion cut
399 c------------------------------------------------------------------------
400 if( (rmax-lmax+1).ge.nclstrp )goto 210 !exits inclusion loop
401
402 if(rstop.eq.0) then !if right cluster morder has not been reached
403 if(value(ir).gt.clinclcut(ir)) then
404 rmax=ir !include a strip on the right
405 else
406 rstop=1 !cluster right end
407 endif
408 endif
409
410 if( (rmax-lmax+1).ge.nclstrp )goto 210 !exits inclusion loop
411
412 if(lstop.eq.0) then !if left cluster morder has not been reached
413 if(value(il).gt.clinclcut(il)) then
414 lmax=il !include a strip on the left
415 else
416 lstop=1 !cluster left end
417 endif
418 endif
419
420 enddo !ends strip inclusion loop
421 goto 211
422 210 continue !jumps here if more than nclstrp have been included
423 c print*,'>>> nclstrp! '
424 211 continue
425 c-----------------------------------------
426 c end of inclusion loop!
427 c-----------------------------------------
428
429 c------------------------------------------------------------------------
430 c adjust the cluster in order to have at least a strip around the seed
431 c------------------------------------------------------------------------
432 if(iseed.eq.lmax.and.lmax.ne.first)then
433 lmax = lmax-1
434 if( (rmax-lmax+1).gt.nclstrp )rmax=rmax-1
435 endif
436 if(iseed.eq.rmax.and.rmax.ne.last )then
437 rmax = rmax+1
438 if( (rmax-lmax+1).gt.nclstrp )lmax=lmax+1
439 endif
440
441 c------------------------------------------------------------------------
442 c adjust the cluster in order to store a minimum number of strips
443 c------------------------------------------------------------------------
444 do while( (rmax-lmax+1).lt.nclstrpmin )
445
446 vl = -99999
447 vr = -99999
448 if(lmax-1.ge.first) vl = value(lmax-1)
449 if(rmax+1.le.last ) vr = value(rmax+1)
450 if(vr.ge.vl) then
451 rmax = rmax+1
452 else
453 lmax = lmax-1
454 endif
455
456 enddo
457
458 c--------------------------------------------------------
459 c store cluster info
460 c--------------------------------------------------------
461 nclstr_view = nclstr_view + 1 !cluster number
462
463 if(nclstr_view.gt.nclstrmax_view) then !too many clusters for the view:
464 c$$$ if(verbose) print*,'Event ',eventn(1),
465 c$$$ $ ': more than ',nclstrmax_view
466 c$$$ $ ,' clusters on view ',iv
467 flag_shower = .true.
468 goto 2000
469 endif
470
471 ladder_view(nclstr_view) = nld(iseed,iv)
472 maxs_view(nclstr_view) = iseed
473 mult_view(nclstr_view) = rmax-lmax+1
474 rmax_view(nclstr_view) = rmax
475 lmax_view(nclstr_view) = lmax
476
477 c$$$ if(rmax-lmax+1.gt.25)
478 c$$$ $ print*,'view ',iv
479 c$$$ $ ,' cl ',nclstr_view,' mult ',rmax-lmax+1
480 c------------------------------------------------------------------------
481 c search for a double peak inside the cluster
482 c------------------------------------------------------------------------
483 inext = rmax+1 !<< index where to start new-cluster search
484
485 vmax = 0
486 vmin = value(iseed)
487 imax = iseed
488 imin = iseed
489 do iss = max(iseed+1,lsat+1),rmax
490 if( value(iss).lt.vmin )then
491 if( imax.ne.iseed )goto 221 !found dowble peek
492 imin = iss
493 vmin = value(iss)
494 else
495 delta = value(iss) - value(imin)
496 cut = sqrt(clinclcut(iss)**2 + clinclcut(imin)**2)
497 if(
498 $ delta.gt.cut .and.
499 $ value(iss).gt.clseedcut(iss).and.
500 $ .true.)then
501 if( value(iss).gt.vmax )then
502 imax = iss
503 vmax = value(iss)
504 else
505 goto 221 !found dowble peek
506 endif
507 endif
508 endif
509 enddo
510 221 continue
511
512 if(imax.gt.iseed)then
513 inext = imax !<< index where to start new-cluster search
514 c$$$ print*,'--- double peek ---'
515 c$$$ print*,(value(ii),ii=lmax,rmax)
516 c$$$ print*,'seed ',iseed,' imin ',imin,' imax ',imax
517 endif
518 c--------------------------------------------------------
519 c
520 c--------------------------------------------------------
521 endif !end possible seed conditio
522 220 continue !jumps here to skip strips left of last seed
523
524 enddo ! end loop on strips
525 enddo !end loop on ladders
526 2000 continue
527 return
528 end
529
530
531 *---***---***---***---***---***---***---***---***
532 *
533 *
534 *
535 *
536 *
537 *---***---***---***---***---***---***---***---***
538
539 subroutine save_cluster(iv)
540 *
541 * (080/2006 Elena Vannuccini)
542 * Save the clusters view by view
543
544 include 'commontracker.f'
545 include 'level1.f'
546 include 'calib.f'
547 include 'common_reduction.f'
548
549 integer CLlength !lunghezza in strip del cluster
550
551 do ic=1,nclstr_view
552
553 nclstr1 = nclstr1+1
554 view(nclstr1) = iv
555 ladder(nclstr1) = ladder_view(ic)
556 maxs(nclstr1) = maxs_view(ic)
557 mult(nclstr1) = mult_view(ic)
558
559 c posizione dell'inizio del cluster nell' array clsignal
560 indstart(nclstr1) = ind
561 c posizione del cluster seed nell'array clsignal
562 indmax(nclstr1) = indstart(nclstr1)
563 $ +( maxs_view(ic) - lmax_view(ic) )
564
565 CLlength = rmax_view(ic) - lmax_view(ic) + 1 !numero di strip salvate
566 totCLlength = totCLlength + CLlength
567 sgnl(nclstr1) = 0
568 do j=lmax_view(ic),rmax_view(ic) !stores sequentially cluter strip values in
569
570 clsignal(ind) = value(j) ! clsignal array
571
572 ivk=nvk(j)
573 ist=nst(j)
574
575 clsigma(ind) = sigma(iv,ivk,ist)
576 cladc(ind) = adc(iv,ivk,ist)
577 clbad(ind) = bad(iv,ivk,ist)
578 c clped(ind) = pedestal(iv,ivk,ist)
579
580 ind=ind+1
581 c if(value(j).gt.0)
582 if(value(j).gt.clinclcut(j))
583 $ sgnl(nclstr1) = sgnl(nclstr1) + value(j) !cluster charge
584 enddo
585
586 c print*,'view ',iv,' -- save_cluster -- nclstr1: '
587 c $ ,nclstr1,maxs(nclstr1),mult(nclstr1),sgnl(nclstr1)
588
589 enddo
590
591 return
592 end
593 *---***---***---***---***---***---***---***---***
594 *
595 *
596 *
597 *
598 *
599 *---***---***---***---***---***---***---***---***
600
601
602 c$$$ subroutine stripmask
603 c$$$
604 c$$$* this routine set va1 and single-strip masks,
605 c$$$* on the basis of the VA1 mask saved in the DB
606 c$$$*
607 c$$$* mask(nviews,nva1_view,nstrips_va1) !strip mask
608 c$$$* mask_vk(nviews,nva1_view) !VA1 mask
609 c$$$*
610 c$$$ include 'commontracker.f'
611 c$$$ include 'level1.f'
612 c$$$ include 'common_reduction.f'
613 c$$$ include 'calib.f'
614 c$$$
615 c$$$* init mask
616 c$$$ do iv=1,nviews
617 c$$$ do ivk=1,nva1_view
618 c$$$ do is=1,nstrips_va1
619 c$$$c mask(iv,ivk,is) = mask_vk(iv,ivk)
620 c$$$ if( mask_vk(iv,ivk) .ne. -1)then
621 c$$$ mask(iv,ivk,is) = 1
622 c$$$ $ * mask_vk(iv,ivk) !from DB
623 c$$$ $ * mask_vk_ev(iv,ivk) !from <SIG>
624 c$$$ $ * mask_vk_run(iv,ivk) !from CN
625 c$$$ else
626 c$$$ mask(iv,ivk,is) = -1
627 c$$$ $ * mask_vk(iv,ivk) !from DB
628 c$$$ $ * mask_vk_ev(iv,ivk) !from CN
629 c$$$ endif
630 c$$$ enddo
631 c$$$ enddo
632 c$$$ enddo
633 c$$$
634 c$$$
635 c$$$ return
636 c$$$ end
637
638 subroutine stripmask(iv,ivk)
639
640 * -----------------------------------------------
641 * this routine set va1 and single-strip masks,
642 * on the basis of the VA1 mask saved in the DB
643 *
644 * mask(nviews,nva1_view,nstrips_va1) !strip mask
645 * mask_vk(nviews,nva1_view) !VA1 mask
646 * -----------------------------------------------
647 include 'commontracker.f'
648 include 'level1.f'
649 include 'common_reduction.f'
650 include 'calib.f'
651
652 * init mask
653 do is=1,nstrips_va1
654 * --------------------------------------------------------
655 * if VA1-mask from DB is 0 or 1, three masks are combined:
656 * - from DB (a-priori mask)
657 * - run-based (chip declared bad on the basis of <SIG>)
658 * - event-based (failure in CN computation)
659 * --------------------------------------------------------
660 if( mask_vk(iv,ivk) .ne. -1)then
661 mask(iv,ivk,is) = 1
662 $ * mask_vk(iv,ivk) !from DB
663 $ * mask_vk_ev(iv,ivk) !from <SIG>
664 $ * mask_vk_run(iv,ivk) !from CN
665 * -----------------------------------------------------------
666 * if VA1-mask from DB is -1 only event-based mask is applied
667 * -----------------------------------------------------------
668 else
669 mask(iv,ivk,is) = -1
670 $ * mask_vk(iv,ivk) !from DB
671 $ * mask_vk_ev(iv,ivk) !from CN
672 endif
673 enddo
674
675
676 return
677 end

  ViewVC Help
Powered by ViewVC 1.1.23