/[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.17 - (show annotations) (download)
Thu Mar 15 12:17:10 2007 UTC (17 years, 9 months ago) by pam-fi
Branch: MAIN
Changes since 1.16: +21 -5 lines
workaround to retrieve clusters + other minor adjustments

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

  ViewVC Help
Powered by ViewVC 1.1.23