/[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.11 - (show annotations) (download)
Fri Oct 27 16:08:19 2006 UTC (18 years, 2 months ago) by pam-fi
Branch: MAIN
Changes since 1.10: +1 -1 lines
new version of the Hough transform

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 * -------------------------------------------------------
27 * STRIP MASK
28 * -------------------------------------------------------
29
30 c call stripmask !called later, after CN computation
31 call init_level1
32
33 c good1 = good0
34 c--------------------------------------------------
35 c check the LEVEL0 event status for missing
36 c sections or DSP alarms
37 c ==> set the variable GOOD1(12)
38 c--------------------------------------------------
39 do iv=1,nviews
40 if(DSPnumber(iv).gt.0.and.DSPnumber(iv).le.12)then
41 c ------------------------
42 c GOOD
43 c ------------------------
44 GOOD1(DSPnumber(iv))=0 !OK
45 c ------------------------
46 c CRC error
47 c ------------------------
48 if(crc(iv).eq.1) then
49 GOOD1(DSPnumber(iv)) = 2
50 goto 18 !next view
51 endif
52 c ------------------------
53 c online-software alarm
54 c ------------------------
55 if(
56 $ fl1(iv).ne.0.or.
57 $ fl2(iv).ne.0.or.
58 $ fl3(iv).ne.0.or.
59 $ fl4(iv).ne.0.or.
60 $ fl5(iv).ne.0.or.
61 $ fl6(iv).ne.0.or.
62 $ fc(iv).ne.0.or.
63 $ DATAlength(iv).eq.0.or.
64 $ .false.)then
65 GOOD1(DSPnumber(iv))=3
66 goto 18
67 endif
68 c ------------------------
69 c DSP-counter jump
70 c ------------------------
71 if(
72 $ eventn_old(iv).ne.0.and. !first event in this file
73 $ eventn(iv).ne.1.and. !first event in run
74 $ good_old(DSPnumber(iv)).ne.0.and. !previous event corrupted
75 $ .true.)then
76
77 if(eventn(iv).ne.(eventn_old(iv)+1))then
78 GOOD1(DSPnumber(iv))=4
79 goto 18
80 endif
81
82 endif
83 c ------------------------
84 18 continue
85 endif
86 enddo
87
88 ngood = 0
89 do iv = 1,nviews
90 eventn_old(iv) = eventn(iv)
91 good_old(iv) = good1(iv)
92 ngood = ngood + good1(iv)
93 enddo
94 c if(ngood.ne.0)print*,'* WARNING * LEVEL0 event status: '
95 c $ ,(good1(i),i=1,nviews)
96 c--------------------------------------------------
97 c read the variable DATATRACKER from LEVEL0
98 c and fill the variable ADC (invertin view 11)
99 c--------------------------------------------------
100 call filladc(iflag)
101 if(iflag.ne.0)then
102 ierror = 220
103 endif
104
105 c--------------------------------------------------
106 c computes common noise for each VA1
107 c (excluding strips with signal,
108 c tagged with the flag CLSTR)
109 c--------------------------------------------------
110 do iv=1,nviews
111 ima=0
112 do ik=1,nva1_view
113 cn(iv,ik) = 0
114 cnrms(iv,ik) = 0
115 cnn(iv,ik) = -1
116 mask_vk_ev(iv,ik)=1
117 iflag=0
118 if(mask_vk(iv,ik).eq.1)call cncomp(iv,ik,iflag)
119 if(iflag.ne.0)then
120 ima=ima+1
121 mask_vk_ev(iv,ik)=0
122 ierror = 220
123 endif
124 enddo
125 100 format(' * WARNING * Event ',i7,' view',i3,': VK MASK ',24i1)
126 if(ima.ne.0.and.debug)write(*,100)eventn(1),iv
127 $ ,(mask_vk_ev(iv,ik),ik=1,nva1_view)
128 enddo
129
130 call stripmask !compute mask(i,j,k), combining mask_vk_ev and mask_vk
131
132 c---------------------------------------------
133 c loops on views, VA1 and strips,
134 c and computes strips signals using
135 c badstrip, pedestals, and
136 c sigma informations from histograms
137 c---------------------------------------------
138 ind=1 !clsignal array index
139
140 do iv=1,nviews !loop on views
141 do is=1,nstrips_view !loop on strips (1)
142 if(mod(iv,2).eq.1) then
143 C=== > Y view
144 value(is)= -(DBLE(adc(iv,nvk(is),nst(is)))
145 $ -cn(iv,nvk(is))-pedestal(iv,nvk(is),nst(is)))
146 $ *mask(iv,nvk(is),nst(is))
147 clseedcut(is)=clcuty*sigma(iv,nvk(is),nst(is))
148 $ *mask(iv,nvk(is),nst(is))
149 clinclcut(is)=incuty*sigma(iv,nvk(is),nst(is))
150 $ *mask(iv,nvk(is),nst(is))
151 sat(is)=0
152 if( adc(iv,nvk(is),nst(is)).lt.adc_saty )sat(is)=1
153 else
154 C=== > X 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)=clcutx*sigma(iv,nvk(is),nst(is))
159 $ *mask(iv,nvk(is),nst(is))
160 clinclcut(is)=incutx*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)).gt.adc_satx )sat(is)=1
164 endif
165 enddo !end loop on strips (1)
166 call search_cluster(iv)
167
168 if(.not.flag_shower)then
169 call save_cluster(iv)
170 else
171 fshower(iv) = 1
172 GOOD1(DSPnumber(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 if(ngood.ne.0)print*,'* WARNING * Event ',eventn(1)
193 $ ,':LEVEL1 event status: '
194 $ ,(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 subroutine stripmask
582
583 * this routine set va1 and single-strip masks,
584 * on the basis of the VA1 mask saved in the DB
585 *
586 * mask(nviews,nva1_view,nstrips_va1) !strip mask
587 * mask_vk(nviews,nva1_view) !VA1 mask
588 *
589 include 'commontracker.f'
590 include 'level1.f'
591 include 'common_reduction.f'
592 include 'calib.f'
593
594 * init mask
595 do iv=1,nviews
596 do ivk=1,nva1_view
597 do is=1,nstrips_va1
598 c mask(iv,ivk,is) = mask_vk(iv,ivk)
599 mask(iv,ivk,is) = mask_vk(iv,ivk) * mask_vk_ev(iv,ivk)
600 enddo
601 enddo
602 enddo
603
604
605 return
606 end
607

  ViewVC Help
Powered by ViewVC 1.1.23