/[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.14 - (show annotations) (download)
Wed Nov 29 10:09:38 2006 UTC (18 years, 1 month ago) by pam-fi
Branch: MAIN
CVS Tags: v2r01
Changes since 1.13: +2 -1 lines
fixed bug on GOOD1 assignment

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

  ViewVC Help
Powered by ViewVC 1.1.23