/[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.19 - (show annotations) (download)
Mon May 14 11:03:06 2007 UTC (17 years, 7 months ago) by pam-fi
Branch: MAIN
CVS Tags: v3r05, v3r06
Changes since 1.18: +4 -0 lines
implemented method to reprocess a track, starting from cluster positions

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

  ViewVC Help
Powered by ViewVC 1.1.23