/[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.13 - (show annotations) (download)
Thu Nov 23 18:51:45 2006 UTC (18 years, 1 month ago) by pam-fi
Branch: MAIN
Changes since 1.12: +57 -19 lines
implemented VA1 mask based on <SIG>

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

  ViewVC Help
Powered by ViewVC 1.1.23