/[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.2 - (hide annotations) (download)
Tue May 30 16:30:37 2006 UTC (18 years, 6 months ago) by pam-fi
Branch: MAIN
CVS Tags: v0r02, v1r00
Changes since 1.1: +19 -178 lines
Error handling from F77 routine / Fixed some bugs with default calibration

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

  ViewVC Help
Powered by ViewVC 1.1.23