/[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.4 - (hide annotations) (download)
Fri Aug 4 08:18:06 2006 UTC (18 years, 4 months ago) by pam-fi
Branch: MAIN
Changes since 1.3: +27 -15 lines
some memory leak bugs fixed + CN computation modified

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.2 integer ierror
22     ierror = 0
23 mocchiut 1.1
24     * -------------------------------------------------------
25     * STRIP MASK
26     * -------------------------------------------------------
27    
28 pam-fi 1.4 c call stripmask !called later, after CN computation
29 mocchiut 1.1 call init_level1
30    
31 pam-fi 1.4 good1 = good0
32 mocchiut 1.1 c--------------------------------------------------
33     c read the variable DATATRACKER from LEVEL0
34     c and fill the variable ADC (inverting view 11)
35     c--------------------------------------------------
36     call filladc(iflag)
37     if(iflag.ne.0)then
38     good1=0
39 pam-fi 1.2 c if(DEBUG)print*,'event ',eventn(1),' >>>>> decode ERROR'
40 pam-fi 1.4 ierror = 220
41 mocchiut 1.1 goto 200
42     endif
43    
44     c--------------------------------------------------
45     c computes common noise for each VA1
46     c (excluding strips affected by signal,
47     c tagged with the flag CLSTR)
48     c--------------------------------------------------
49     do iv=1,nviews
50     do ik=1,nva1_view
51 pam-fi 1.4 cn(iv,ik)=0
52     mask_vk_ev(iv,ik)=1
53 pam-fi 1.2 iflag=0
54     if(mask_vk(iv,ik).eq.1)call cncomp(iv,ik,iflag)
55 pam-fi 1.4 c if(iflag.ne.0)good1=0
56     if(iflag.ne.0)then
57     mask_vk_ev(iv,ik)=0
58     ierror = 220
59     endif
60 mocchiut 1.1 enddo
61     enddo
62 pam-fi 1.4 c if(good1.eq.0)then
63     c ierror = 220
64     c endif
65 mocchiut 1.1
66 pam-fi 1.4 call stripmask !compute mask(i,j,k)
67 mocchiut 1.1 c---------------------------------------------
68     c loops on views, VA1 and strips,
69     c and computes strips signals using
70     c badstrip, pedestals, and
71     c sigma informations from histograms
72     c---------------------------------------------
73     flag_shower = .false.
74     ind=1 !clsignal array index
75     do iv=1,nviews !loop on views
76     do is=1,nstrips_view !loop on strips (1)
77     if(mod(iv,2).eq.1) then
78     C=== > Y view
79     value(is)= -(DBLE(adc(iv,nvk(is),nst(is)))
80     $ -cn(iv,nvk(is))-pedestal(iv,nvk(is),nst(is)))
81     $ *mask(iv,nvk(is),nst(is))
82     clseedcut(is)=clcuty*sigma(iv,nvk(is),nst(is))
83     $ *mask(iv,nvk(is),nst(is))
84     clinclcut(is)=incuty*sigma(iv,nvk(is),nst(is))
85     $ *mask(iv,nvk(is),nst(is))
86     ccc print*,"value(",is,")(reduction)= ",value(is)
87     else
88     C=== > X view
89     value(is)= (DBLE(adc(iv,nvk(is),nst(is)))
90     $ -cn(iv,nvk(is))-pedestal(iv,nvk(is),nst(is)))
91     $ *mask(iv,nvk(is),nst(is))
92     clseedcut(is)=clcutx*sigma(iv,nvk(is),nst(is))
93     $ *mask(iv,nvk(is),nst(is))
94     clinclcut(is)=incutx*sigma(iv,nvk(is),nst(is))
95     $ *mask(iv,nvk(is),nst(is))
96     endif
97 pam-fi 1.4 c$$$ print*,iv,is,' --- ',adc(iv,nvk(is),nst(is)),cn(iv,nvk(is))
98     c$$$ $ ,pedestal(iv,nvk(is),nst(is)),value(is)
99     c$$$ $ ,sigma(iv,nvk(is),nst(is))
100     c if(value(is).gt.clseedcut(is))
101     c $ print*,iv,is,' --- (ADC_PED_CN) ',value(is),clseedcut(is)
102 mocchiut 1.1 enddo !end loop on strips (1)
103     call search_cluster(iv)
104     if(flag_shower.eqv..true.)then
105     call init_level1
106     good1=0
107     goto 200 !jump to next event
108     endif
109     enddo ! end loop on views
110     do iv=1,nviews
111     do ik=1,nva1_view
112     cnev(iv,ik)=cn(iv,ik) !assigns computed CN to ntuple variables
113 pam-fi 1.4 cnevflag(iv,ik)=cnflag(iv,ik) !assigns computed CN to ntuple variables
114 mocchiut 1.1 ccc print*,"cnev(",iv,",",ik,")(reduction)= ",cnev(iv,ik)
115     enddo
116     enddo
117     C---------------------------------------------
118     C come here if GOOD1=0
119     C or the event has too many clusters
120     C---------------------------------------------
121     200 continue
122     c------------------------------------------------------------------------
123 pam-fi 1.2 c
124 mocchiut 1.1 c closes files and exits
125 pam-fi 1.2 c
126 mocchiut 1.1 c------------------------------------------------------------------------
127 pam-fi 1.2 RETURN
128     END
129 mocchiut 1.1
130     ***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***
131     *
132     *
133     *
134     *
135     *
136     *
137     *
138     *
139     *
140     ***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***
141    
142    
143     subroutine init_level1
144    
145     include 'commontracker.f'
146     include 'level1.f'
147     include 'level0.f'
148    
149     good1=0
150     nclstr1=0
151     totCLlength=0
152     do ic=1,nclstrmax
153     view(ic)=0
154     ladder(ic)=0
155     indstart(ic)=0
156     indmax(ic)=0
157     maxs(ic)=0
158     mult(ic)=0
159     dedx(ic)=0
160     enddo
161     do id=1,maxlength !???
162     clsignal(id)=0.
163     enddo
164     do iv=1,nviews
165     c crc1(iv)=0
166     do ik=1,nva1_view
167     cnev(iv,ik)=0
168     enddo
169     enddo
170    
171     return
172     end
173     *---***---***---***---***---***---***---***---***
174     *
175     *
176     *
177     *
178     *
179     *---***---***---***---***---***---***---***---***
180    
181     subroutine search_cluster(iv)
182    
183     include 'commontracker.f'
184     include 'common_reduction.f'
185     include 'level0.f'
186     include 'level1.f'
187     include 'calib.f'
188    
189    
190    
191     c local variables
192     integer rmax,lmax !estremi del cluster
193     integer rstop,lstop !per decidere quali strip includere nel cluster
194     ! oltre il seed
195     integer first,last,diff !per includere le strip giuste... !???
196    
197     integer multtemp !temporary multiplicity variable
198    
199     integer CLlength !lunghezza in strip del cluster
200    
201     external nst
202    
203     c------------------------------------------------------------------------
204     c looks for clusters on each view
205     C : CERCO STRIP SOPRA CLSEEDCUT, POI SCORRO A DX FINCHE'
206     c NON TROVO
207     C STRIP PIU' BASSA (in segnale/rumore)
208     C => L'ULTIMA DELLA SERIE CRESCENTE
209     C (LA PIU' ALTA) E' IL
210     C CLUSTER SEED. POI SCORRO A SX E DX INCLUDENDO TUTTE
211     C LE STRIP (FINO A 17 AL
212     C MAX) CHE SUPERANO CLINCLCUT.
213     C QUANDO CERCO IL CLUSTER SEED SUCCESSIVO SALTO LA STRIP
214     C ADIACENTE A DESTRA
215     C DELL'ULTIMO CLUSTER SEED (CHE SARA' NECESSARIAMENTE
216     C PIU' BASSA) E PRENDO
217     C COME SEED UNA STRIP SOLO SE IL SUO SEGNALE E'
218     C MAGGIORE DI QUELLO DELLA STRIP
219     C PRECEDENTE (PRATICAMENTE PER EVITARE CHE L'ULTIMA
220     C STRIP DI UN GRUPPO DI STRIP
221     C TUTTE SOPRA IL CLSEEDCUT VENGA AUTOMATICAMENTE PRESA
222     C COME SEED... DEVE ESSERE
223     C PRESA SOLO SE IL CLUSTER E' DOUBLE PEAKED...)
224     c------------------------------------------------------------------------
225     c 6 ottobre 2003
226     c Elena: CLSEEDCUT = 7 (old value 10)
227     c Elena: CLINCLCUT = 4 (old value 5)
228    
229     iseed=-999 !cluster seed index initialization
230    
231     do jl=1,nladders_view !1..3 !loops on ladders
232     first=1+nstrips_ladder*(jl-1) !1,1025,2049
233     last=nstrips_ladder*jl !1024,2048,3072
234     c X views have 1018 strips instead of 1024
235     if(mod(iv,2).eq.0) then
236     first=first+3
237     last=last-3
238     endif
239     do is=first,last !loop on strips in each ladder
240     if(is.le.iseed+1) goto 220
241     c-----------------------------------------
242     c after a cluster seed as been found,
243     c look for next one skipping one strip on the right
244     c (i.e. look for double peak cluster)
245     c-----------------------------------------
246     if(is.ne.first) then
247     if(value(is).le.value(is-1)) goto 220
248     endif
249     c-----------------------------------------
250     c skips cluster seed
251     c finding if strips values are descreasing (a strip
252     c can be a cluster seed only if previous strip value
253     c is lower)
254     c-----------------------------------------
255     if(value(is).gt.clseedcut(is)) then
256     ccc print*,"value(",is,")=",value(is),
257     ccc $ " .gt.clseedcut(",is,")=",clseedcut(is)
258     c-----------------------------------------
259     c possible SEED...
260     c-----------------------------------------
261     itemp=is
262     if(itemp.eq.last) goto 230 !estremo...
263     do while(value(itemp)
264     $ /sigma(iv,nvk(itemp),nst(itemp))
265     $ .le.value(itemp+1)
266     $ /sigma(iv,nvk(itemp+1),nst(itemp+1))) !BIAS: aggiustare il caso uguale!???
267     itemp=itemp+1
268     if(itemp.eq.last) goto 230 !stops if reaches last strip
269     enddo ! of the ladder
270     230 continue
271     c-----------------------------------------
272     c fownd SEED!!!
273     c-----------------------------------------
274     iseed=itemp
275     c----------------------------------------------------------
276     c after finding a cluster seed, checks also adjacent strips,
277     C and marks the ones exceeding clinclcut
278     c----------------------------------------------------------
279     ir=iseed !indici destro
280     il=iseed ! e sinistro
281    
282     rmax=ir !estremo destro del cluster
283     lmax=il ! e sinistro
284    
285     rstop=0 !initialize flags used to exit from
286     lstop=0 ! inclusion loop
287    
288     do while(lstop.eq.0.or.rstop.eq.0) !shifts left and right from
289     ir=ir+1 !position index for strips on right side of
290     ! cluster seed
291     il=il-1 !and for left side
292     c------------------------------------------------------------------------
293     c checks for last or first strip of the ladder
294     c------------------------------------------------------------------------
295     if(ir.gt.last) then !when index goes beyond last strip
296     rstop=1 ! of the ladder, change rstop flag in order
297     ! to "help" exiting from loop
298     endif
299    
300     if(il.lt.first) then !idem when index goes beyond
301     lstop=1 ! first strip of the ladder
302     endif
303    
304     c------------------------------------------------------------------------
305     c check for clusters including more than nclstrp strips
306     c------------------------------------------------------------------------
307     if((rmax-lmax+1).ge.nclstrp) then
308     goto 210 !exits inclusion loop:
309     ! lmax and rmax maintain last value
310     ! NB .ge.!???
311     endif
312     c------------------------------------------------------------------------
313     c marks strips exceeding inclusion cut
314     c------------------------------------------------------------------------
315     if(rstop.eq.0) then !if last strip of the ladder or last
316     ! over-cut strip has not been reached
317     if(value(ir).gt.clinclcut(ir)) then !puts in rmax the
318     rmax=ir ! last right over-cut strip
319     else
320     rstop=1 !otherwise cluster ends on right and rstop
321     endif ! flag=1 signals it
322     endif
323     if(lstop.eq.0) then
324     if(value(il).gt.clinclcut(il)) then
325     lmax=il
326     else
327     lstop=1
328     endif
329     endif
330    
331     enddo !ends strip inclusion loop
332     210 continue !jumps here if more than nclstrp have been included
333    
334     multtemp=rmax-lmax+1 !stores multiplicity in temp
335     ! variable. NB rmax and lmax can change later in
336     ! order to include enough strips to calculate eta3
337     ! and eta4. so mult is not always equal to cllength
338     c------------------------------------------------------------------------
339     c NB per essere sicuro di poter calcolare eta3 e eta4 devo includere
340     c sempre e comunque le 2 strip adiacenti al cluster seed e quella
341     c adiacente ulteriore dalla parte della piu' alta fra queste due
342     c (vedi oltre...)!???
343     c------------------------------------------------------------------------
344    
345     c nel caso di estremi del ladder...!???
346    
347     c ho meno di 4 strip nel cluster --> se sono sui bordi o quasi del ladder
348     c costruisco il cluster ad hoc e poi esco, se non sono sui bordi o quasi
349     c vado oltre (aggiungero' quindi strip a sx e dx in modo da poter calcolare
350     c eta3e4)
351     if((rmax-lmax+1).lt.4) then
352    
353     if(iseed.eq.first) then !estremi...
354     rmax=iseed+2 !NB in questo modo puo' anche capitare di
355     lmax=iseed ! includere strip sotto taglio di inclusione
356     goto 250 ! che non serviranno per eta3e4!???
357     endif
358    
359     if(iseed.eq.last) then !estremi...
360     rmax=iseed
361     lmax=iseed-2 !NB 2 e non 3, perche' altrimenti sarei in
362     goto 250 ! ((rmax-lmax+1).lt.4).eq.false. !???
363     endif !NMB questo e' l'unico caso di cllength=3!???
364    
365     if(iseed.eq.first+1) then !quasi estremi...
366     rmax=iseed+2
367     lmax=iseed-1
368     goto 250
369     endif
370     if(iseed.eq.last-1) then
371     rmax=iseed+1
372     lmax=iseed-2
373     goto 250
374     endif
375     c se ho 4 o piu' strip --> se sono sui bordi esco, se sono sui quasi bordi
376     c includo la strip del bordo
377     else
378    
379     if(iseed.eq.first) goto 250 !estremi... non includo altro
380     if(iseed.eq.last) goto 250
381     if(iseed.eq.first+1) then !quasi estremi... mi assicuro di
382     lmax=first ! avere le strip adiacenti al seed
383     if((rmax-lmax+1).gt.nclstrp) rmax=rmax-1 !NB effetto
384     goto 250 ! coperta: se la lunghezza del cluster era gia'
385     endif ! al limite (nclstrp), per poter aggiungere questa
386     ! strip a sinistra devo toglierne una a destra...!???
387     if(iseed.eq.last-1) then
388     rmax=last
389     if((rmax-lmax+1).gt.nclstrp) lmax=lmax+1
390     goto 250
391     endif
392     endif
393     c------------------------------------------------------------------------
394     c be sure to include in the cluster the cluster seed with its 2 adjacent
395     c strips, and the one adjacent to the greatest between this two strip, as the
396     c fourth one. if the strips have the same value (!) the fourth one is chosen
397     c as the one having the greatest value between the second neighbors
398     c------------------------------------------------------------------------
399     if(value(iseed+1).eq.value(iseed-1)) then
400     if(value(iseed+2).ge.value(iseed-2)) then !??? qui cmq c'e'
401     diff=(iseed+2)-rmax
402     if(diff.gt.0) then
403     rmax=rmax+diff
404     if((rmax-lmax+1).gt.nclstrp) then
405     lmax=rmax-nclstrp+1
406     endif
407     endif
408     diff=(iseed-1)-lmax
409     if(diff.lt.0) then
410     lmax=lmax+diff
411     if((rmax-lmax+1).gt.nclstrp) then
412     rmax=lmax+nclstrp-1
413     endif
414     endif
415     else
416     diff=(iseed-2)-lmax
417     if(diff.lt.0) then
418     lmax=lmax+diff
419     if((rmax-lmax+1).gt.nclstrp) then
420     rmax=lmax+nclstrp-1
421     endif
422     endif
423     diff=(iseed+1)-rmax
424     if(diff.gt.0) then
425     rmax=rmax+diff
426     if((rmax-lmax+1).gt.nclstrp) then
427     lmax=rmax-nclstrp+1
428     endif
429     endif
430     endif
431     elseif(value(iseed+1).gt.value(iseed-1)) then
432     c !??? sposto il limite del cluster a destra per includere sempre le strip
433     c necessarie al calcolo di eta-i
434     c se il cluster diventa troppo lungo lo accorcio a sinistra per avere non piu'
435     c di nclstrp (in questo caso sono sicuro di aver gia' incluso le strip
436     c necessarie al calcolo di eta-i a sinistra, quindi se voglio posso uscire)
437     diff=(iseed+2)-rmax
438     if(diff.gt.0) then
439     rmax=rmax+diff
440     if((rmax-lmax+1).gt.nclstrp) then
441     lmax=rmax-nclstrp+1
442     c goto 250
443     endif
444     endif
445     diff=(iseed-1)-lmax
446     if(diff.lt.0) then
447     lmax=lmax+diff
448     if((rmax-lmax+1).gt.nclstrp) then
449     rmax=lmax+nclstrp-1
450     c goto 250 !inutile!???
451     endif
452     endif
453     else
454     diff=(iseed-2)-lmax
455     if(diff.lt.0) then
456     lmax=lmax+diff
457     if((rmax-lmax+1).gt.nclstrp) then
458     rmax=lmax+nclstrp-1
459     c goto 250
460     endif
461     endif
462     diff=(iseed+1)-rmax
463     if(diff.gt.0) then
464     rmax=rmax+diff
465     if((rmax-lmax+1).gt.nclstrp) then
466     lmax=rmax-nclstrp+1
467     c goto 250 !inutile!???
468     endif
469     endif
470     endif
471     250 continue
472    
473     c--------------------------------------------------------
474 pam-fi 1.4 c fills cluster variables
475 mocchiut 1.1 c--------------------------------------------------------
476     nclstr1=nclstr1+1 !cluster number
477     ccc print*,nclstr1,multtemp
478     if(nclstr1.gt.nclstrmax) then !too many clusters for the event:
479 pam-fi 1.4 if(verbose)print*,'Event ',eventn(1),
480     $ ': more than ',nclstrmax,' clusters'
481 mocchiut 1.1 good1=0 ! event
482     nclstr1=0
483     totCLlength=0
484     flag_shower = .true.
485     goto 2000
486     endif
487     view(nclstr1)=iv !vista del cluster
488     ladder(nclstr1)=nld(iseed,iv) !ladder a cui appartiene il cluster seed
489     maxs(nclstr1)=iseed !strip del cluster seed
490     mult(nclstr1)=multtemp !molteplicita'
491    
492     indstart(nclstr1)=ind !posizione dell'inizio del cluster nell'
493     ! array clsignal
494     indmax(nclstr1)=indstart(nclstr1)+(iseed-lmax) !posizione del
495     ! cluster seed nell'array clsignal
496    
497     CLlength=rmax-lmax+1 !numero di strip del cluster
498     totCLlength=totCLlength+CLlength
499     dedx(nclstr1)=0
500     do j=lmax,rmax !stores sequentially cluter strip values in
501     clsignal(ind)=value(j) ! clsignal array
502     ind=ind+1
503     c if(value(j).gt.0)
504     if(value(j).gt.clinclcut(j))
505     $ dedx(nclstr1)=dedx(nclstr1)+value(j) !cluster charge
506     enddo
507     c--------------------------------------------------------
508 pam-fi 1.2 c
509 mocchiut 1.1 c--------------------------------------------------------
510     endif !end possible seed conditio
511     220 continue !jumps here to skip strips left of last seed
512    
513     enddo ! end loop on strips
514     enddo !end loop on ladders
515     2000 continue
516     return
517     end
518    
519    
520     *---***---***---***---***---***---***---***---***
521     *
522     *
523     *
524     *
525     *
526     *---***---***---***---***---***---***---***---***
527    
528    
529     subroutine stripmask
530    
531     * this routine set va1 and single-strip masks,
532     * on the basis of the VA1 mask saved in the DB
533     *
534     * mask(nviews,nva1_view,nstrips_va1) !strip mask
535     * mask_vk(nviews,nva1_view) !VA1 mask
536     *
537     include 'commontracker.f'
538 pam-fi 1.4 c include 'level1.f'
539     include 'common_reduction.f'
540 mocchiut 1.1 include 'calib.f'
541    
542     * init mask
543     do iv=1,nviews
544     do ivk=1,nva1_view
545     do is=1,nstrips_va1
546 pam-fi 1.4 c mask(iv,ivk,is) = mask_vk(iv,ivk)
547     mask(iv,ivk,is) = mask_vk(iv,ivk) * mask_vk_ev(iv,ivk)
548 mocchiut 1.1 enddo
549     enddo
550     enddo
551    
552    
553     return
554     end
555    

  ViewVC Help
Powered by ViewVC 1.1.23