/[PAMELA software]/DarthVader/TrackerLevel2/src/F77/reductionflight.f
ViewVC logotype

Annotation of /DarthVader/TrackerLevel2/src/F77/reductionflight.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.10 - (hide annotations) (download)
Thu Oct 26 16:22:38 2006 UTC (18 years, 1 month ago) by pam-fi
Branch: MAIN
Changes since 1.9: +167 -330 lines
fitting methods

1 mocchiut 1.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 pam-fi 1.2 subroutine reductionflight(ierror)
14 mocchiut 1.1
15     include 'commontracker.f'
16     include 'level0.f'
17     include 'level1.f'
18     include 'common_reduction.f'
19     include 'calib.f'
20    
21 pam-fi 1.6 data eventn_old/nviews*0/
22    
23 pam-fi 1.2 integer ierror
24     ierror = 0
25 mocchiut 1.1
26     * -------------------------------------------------------
27     * STRIP MASK
28     * -------------------------------------------------------
29    
30 pam-fi 1.4 c call stripmask !called later, after CN computation
31 mocchiut 1.1 call init_level1
32    
33 pam-fi 1.6 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 pam-fi 1.9 GOOD1(DSPnumber(iv))=0 !OK
45 pam-fi 1.6 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 mocchiut 1.1 c--------------------------------------------------
97     c read the variable DATATRACKER from LEVEL0
98 pam-fi 1.6 c and fill the variable ADC (invertin view 11)
99 mocchiut 1.1 c--------------------------------------------------
100     call filladc(iflag)
101     if(iflag.ne.0)then
102 pam-fi 1.6 ierror = 220
103 mocchiut 1.1 endif
104    
105     c--------------------------------------------------
106     c computes common noise for each VA1
107 pam-fi 1.10 c (excluding strips with signal,
108 mocchiut 1.1 c tagged with the flag CLSTR)
109     c--------------------------------------------------
110     do iv=1,nviews
111 pam-fi 1.8 ima=0
112     do ik=1,nva1_view
113     cn(iv,ik) = 0
114 pam-fi 1.10 cnrms(iv,ik) = 0
115 pam-fi 1.8 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 pam-fi 1.9 100 format(' * WARNING * Event ',i7,' view',i3,': VK MASK ',24i1)
126 pam-fi 1.10 if(ima.ne.0.and.debug)write(*,100)eventn(1),iv
127 pam-fi 1.8 $ ,(mask_vk_ev(iv,ik),ik=1,nva1_view)
128 mocchiut 1.1 enddo
129    
130 pam-fi 1.5 call stripmask !compute mask(i,j,k), combining mask_vk_ev and mask_vk
131 pam-fi 1.10
132 mocchiut 1.1 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 pam-fi 1.5
140 mocchiut 1.1 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 pam-fi 1.10 sat(is)=0
152     if( adc(iv,nvk(is),nst(is)).lt.adc_saty )sat(is)=1
153 mocchiut 1.1 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 pam-fi 1.10 sat(is)=0
163     if( adc(iv,nvk(is),nst(is)).gt.adc_satx )sat(is)=1
164 mocchiut 1.1 endif
165     enddo !end loop on strips (1)
166     call search_cluster(iv)
167 pam-fi 1.10
168 pam-fi 1.5 if(.not.flag_shower)then
169     call save_cluster(iv)
170     else
171     fshower(iv) = 1
172 pam-fi 1.10 GOOD1(DSPnumber(iv)) = 11
173 mocchiut 1.1 endif
174     enddo ! end loop on views
175     do iv=1,nviews
176     do ik=1,nva1_view
177 pam-fi 1.10 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 mocchiut 1.1 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 pam-fi 1.6
188     ngood = 0
189     do iv = 1,nviews
190     ngood = ngood + good1(iv)
191     enddo
192 pam-fi 1.7 if(ngood.ne.0)print*,'* WARNING * Event ',eventn(1)
193     $ ,':LEVEL1 event status: '
194 pam-fi 1.6 $ ,(good1(i),i=1,nviews)
195 mocchiut 1.1 c------------------------------------------------------------------------
196 pam-fi 1.2 c
197 mocchiut 1.1 c closes files and exits
198 pam-fi 1.2 c
199 mocchiut 1.1 c------------------------------------------------------------------------
200 pam-fi 1.2 RETURN
201     END
202 mocchiut 1.1
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 pam-fi 1.6 c good1 = 0
223     do iv=1,12
224     good1(iv) = 1 !missing packet
225     enddo
226 pam-fi 1.5 nclstr1 = 0
227     totCLlength = 0
228 mocchiut 1.1 do ic=1,nclstrmax
229 pam-fi 1.5 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 mocchiut 1.1 enddo
239     do id=1,maxlength !???
240 pam-fi 1.5 clsignal(id) = 0.
241     clsigma(id) = 0.
242     cladc(id) = 0.
243     clbad(id) = 0.
244 mocchiut 1.1 enddo
245     do iv=1,nviews
246     c crc1(iv)=0
247     do ik=1,nva1_view
248 pam-fi 1.5 cnev(iv,ik) = 0
249     cnnev(iv,ik) = 0
250 mocchiut 1.1 enddo
251 pam-fi 1.5 fshower(iv) = 0
252 mocchiut 1.1 enddo
253    
254     return
255     end
256 pam-fi 1.10
257 mocchiut 1.1 *---***---***---***---***---***---***---***---***
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 pam-fi 1.5 include 'common_reduction.f'
273 mocchiut 1.1
274    
275     c local variables
276     integer rmax,lmax !estremi del cluster
277 pam-fi 1.10 integer rstop,lstop
278     integer first,last
279     integer fsat,lsat
280 mocchiut 1.1
281     external nst
282    
283 pam-fi 1.10 iseed=-999 !cluster seed index initialization
284 mocchiut 1.1
285 pam-fi 1.10 inext=-999 !index where to start new cluster search
286 mocchiut 1.1
287 pam-fi 1.10 flag_shower = .false.
288 pam-fi 1.5 nclstr_view=0
289    
290 pam-fi 1.10 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 mocchiut 1.1 if(mod(iv,2).eq.0) then
297     first=first+3
298     last=last-3
299     endif
300 pam-fi 1.6
301 mocchiut 1.1 do is=first,last !loop on strips in each ladder
302 pam-fi 1.6
303 pam-fi 1.10 c---------------------------------------------
304     c new-cluster search starts at index inext
305     c---------------------------------------------
306     if(is.lt.inext) goto 220 ! next strip
307 pam-fi 1.5
308 mocchiut 1.1 if(value(is).gt.clseedcut(is)) then
309     c-----------------------------------------
310     c possible SEED...
311     c-----------------------------------------
312 pam-fi 1.10 itemp = is
313     fsat = 0 ! first saturated strip
314     lsat = 0 ! last saturated strip
315 mocchiut 1.1 if(itemp.eq.last) goto 230 !estremo...
316 pam-fi 1.10 c ------------------------
317     c search for first maximum
318     c ------------------------
319 pam-fi 1.6 do while(
320     $ value(itemp).le.value(itemp+1)
321     $ .and.value(itemp+1).gt.clseedcut(itemp+1))
322 pam-fi 1.10 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 mocchiut 1.1 enddo ! of the ladder
326     230 continue
327 pam-fi 1.10 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 mocchiut 1.1 c fownd SEED!!!
343 pam-fi 1.10 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 mocchiut 1.1 c after finding a cluster seed, checks also adjacent strips,
354 pam-fi 1.10 C and tags the ones exceeding clinclcut
355     c---------------------------------------------------------------
356 mocchiut 1.1 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 pam-fi 1.10
367    
368     ir=ir+1 !index for right side
369     il=il-1 !index for left side
370 mocchiut 1.1 c------------------------------------------------------------------------
371     c checks for last or first strip of the ladder
372     c------------------------------------------------------------------------
373 pam-fi 1.10 if( ir.gt.last ) rstop = 1
374     if( il.lt.first ) lstop = 1
375 mocchiut 1.1
376     c------------------------------------------------------------------------
377 pam-fi 1.10 c add strips exceeding inclusion cut
378 mocchiut 1.1 c------------------------------------------------------------------------
379 pam-fi 1.10 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 mocchiut 1.1 else
385 pam-fi 1.10 rstop=1 !cluster right end
386     endif
387 mocchiut 1.1 endif
388 pam-fi 1.10
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 mocchiut 1.1 if(value(il).gt.clinclcut(il)) then
393 pam-fi 1.10 lmax=il !include a strip on the left
394 mocchiut 1.1 else
395 pam-fi 1.10 lstop=1 !cluster left end
396 mocchiut 1.1 endif
397     endif
398    
399     enddo !ends strip inclusion loop
400 pam-fi 1.10 goto 211
401 mocchiut 1.1 210 continue !jumps here if more than nclstrp have been included
402 pam-fi 1.10 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 mocchiut 1.1 endif
419 pam-fi 1.10
420 mocchiut 1.1 c------------------------------------------------------------------------
421 pam-fi 1.10 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 mocchiut 1.1 endif
434 pam-fi 1.10
435     enddo
436 mocchiut 1.1
437     c--------------------------------------------------------
438 pam-fi 1.10 c store cluster info
439 mocchiut 1.1 c--------------------------------------------------------
440 pam-fi 1.5 nclstr_view = nclstr_view + 1 !cluster number
441 pam-fi 1.10
442 pam-fi 1.5 if(nclstr_view.gt.nclstrmax_view) then !too many clusters for the view:
443 pam-fi 1.10 c$$$ if(verbose) print*,'Event ',eventn(1),
444     c$$$ $ ': more than ',nclstrmax_view
445     c$$$ $ ,' clusters on view ',iv
446 mocchiut 1.1 flag_shower = .true.
447     goto 2000
448     endif
449 pam-fi 1.5
450 pam-fi 1.10 ladder_view(nclstr_view) = nld(iseed,iv)
451     maxs_view(nclstr_view) = iseed
452     mult_view(nclstr_view) = rmax-lmax+1
453 pam-fi 1.5 rmax_view(nclstr_view) = rmax
454     lmax_view(nclstr_view) = lmax
455    
456 pam-fi 1.10 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 dowble 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 mocchiut 1.1 c--------------------------------------------------------
498 pam-fi 1.2 c
499 mocchiut 1.1 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 pam-fi 1.5 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 mocchiut 1.1
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 pam-fi 1.5 include 'level1.f'
591 pam-fi 1.4 include 'common_reduction.f'
592 mocchiut 1.1 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 pam-fi 1.4 c mask(iv,ivk,is) = mask_vk(iv,ivk)
599     mask(iv,ivk,is) = mask_vk(iv,ivk) * mask_vk_ev(iv,ivk)
600 mocchiut 1.1 enddo
601     enddo
602     enddo
603    
604    
605     return
606     end
607    

  ViewVC Help
Powered by ViewVC 1.1.23