*************************************************************************
*
* Program reductionflight.f
*
* - reads readraw.f output files: LEVEL0 ntuple, and ped, sig and bad histograms
* - decodes raw data (DATATRACKER) using DSP ped, sig and bad values
* - looks for clusters information using ped, sig and bad values from
* DSP histograms
* - fills LEVEL1 ntuple
*
*************************************************************************
subroutine reductionflight(ierror)
include 'commontracker.f'
include 'level0.f'
include 'level1.f'
include 'common_reduction.f'
include 'calib.f'
integer ierror
ierror = 0
* -------------------------------------------------------
* STRIP MASK
* -------------------------------------------------------
call stripmask
call init_level1
good1=good0
c--------------------------------------------------
c read the variable DATATRACKER from LEVEL0
c and fill the variable ADC (inverting view 11)
c--------------------------------------------------
call filladc(iflag)
if(iflag.ne.0)then
good1=0
c if(DEBUG)print*,'event ',eventn(1),' >>>>> decode ERROR'
ierror = 220
goto 200
endif
c--------------------------------------------------
c computes common noise for each VA1
c (excluding strips affected by signal,
c tagged with the flag CLSTR)
c--------------------------------------------------
do iv=1,nviews
do ik=1,nva1_view
cn(iv,ik)=0 !initializes cn variable
iflag=0
if(mask_vk(iv,ik).eq.1)call cncomp(iv,ik,iflag)
if(iflag.ne.0)good1=0
enddo
enddo
if(good1.eq.0)then
ierror = 220
c if(WARNING)
c $ print*,' WARNING - cncomp: CN computation failure '
endif
c---------------------------------------------
c loops on views, VA1 and strips,
c and computes strips signals using
c badstrip, pedestals, and
c sigma informations from histograms
c---------------------------------------------
flag_shower = .false.
ind=1 !clsignal array index
do iv=1,nviews !loop on views
do is=1,nstrips_view !loop on strips (1)
if(mod(iv,2).eq.1) then
C=== > Y view
value(is)= -(DBLE(adc(iv,nvk(is),nst(is)))
$ -cn(iv,nvk(is))-pedestal(iv,nvk(is),nst(is)))
$ *mask(iv,nvk(is),nst(is))
clseedcut(is)=clcuty*sigma(iv,nvk(is),nst(is))
$ *mask(iv,nvk(is),nst(is))
clinclcut(is)=incuty*sigma(iv,nvk(is),nst(is))
$ *mask(iv,nvk(is),nst(is))
ccc print*,"value(",is,")(reduction)= ",value(is)
else
C=== > X view
value(is)= (DBLE(adc(iv,nvk(is),nst(is)))
$ -cn(iv,nvk(is))-pedestal(iv,nvk(is),nst(is)))
$ *mask(iv,nvk(is),nst(is))
clseedcut(is)=clcutx*sigma(iv,nvk(is),nst(is))
$ *mask(iv,nvk(is),nst(is))
clinclcut(is)=incutx*sigma(iv,nvk(is),nst(is))
$ *mask(iv,nvk(is),nst(is))
endif
enddo !end loop on strips (1)
call search_cluster(iv)
if(flag_shower.eqv..true.)then
call init_level1
good1=0
goto 200 !jump to next event
endif
enddo ! end loop on views
do iv=1,nviews
do ik=1,nva1_view
cnev(iv,ik)=cn(iv,ik) !assigns computed CN to ntuple variables
ccc print*,"cnev(",iv,",",ik,")(reduction)= ",cnev(iv,ik)
enddo
enddo
C---------------------------------------------
C come here if GOOD1=0
C or the event has too many clusters
C---------------------------------------------
200 continue
c------------------------------------------------------------------------
c
c closes files and exits
c
c------------------------------------------------------------------------
RETURN
END
***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***
*
*
*
*
*
*
*
*
*
***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***
subroutine init_level1
include 'commontracker.f'
include 'level1.f'
include 'level0.f'
good1=0
nclstr1=0
totCLlength=0
do ic=1,nclstrmax
view(ic)=0
ladder(ic)=0
indstart(ic)=0
indmax(ic)=0
maxs(ic)=0
mult(ic)=0
dedx(ic)=0
enddo
do id=1,maxlength !???
clsignal(id)=0.
enddo
do iv=1,nviews
c crc1(iv)=0
do ik=1,nva1_view
cnev(iv,ik)=0
enddo
enddo
return
end
*---***---***---***---***---***---***---***---***
*
*
*
*
*
*---***---***---***---***---***---***---***---***
subroutine search_cluster(iv)
include 'commontracker.f'
include 'common_reduction.f'
include 'level0.f'
include 'level1.f'
include 'calib.f'
c local variables
integer rmax,lmax !estremi del cluster
integer rstop,lstop !per decidere quali strip includere nel cluster
! oltre il seed
integer first,last,diff !per includere le strip giuste... !???
integer multtemp !temporary multiplicity variable
integer CLlength !lunghezza in strip del cluster
external nst
c------------------------------------------------------------------------
c looks for clusters on each view
C : CERCO STRIP SOPRA CLSEEDCUT, POI SCORRO A DX FINCHE'
c NON TROVO
C STRIP PIU' BASSA (in segnale/rumore)
C => L'ULTIMA DELLA SERIE CRESCENTE
C (LA PIU' ALTA) E' IL
C CLUSTER SEED. POI SCORRO A SX E DX INCLUDENDO TUTTE
C LE STRIP (FINO A 17 AL
C MAX) CHE SUPERANO CLINCLCUT.
C QUANDO CERCO IL CLUSTER SEED SUCCESSIVO SALTO LA STRIP
C ADIACENTE A DESTRA
C DELL'ULTIMO CLUSTER SEED (CHE SARA' NECESSARIAMENTE
C PIU' BASSA) E PRENDO
C COME SEED UNA STRIP SOLO SE IL SUO SEGNALE E'
C MAGGIORE DI QUELLO DELLA STRIP
C PRECEDENTE (PRATICAMENTE PER EVITARE CHE L'ULTIMA
C STRIP DI UN GRUPPO DI STRIP
C TUTTE SOPRA IL CLSEEDCUT VENGA AUTOMATICAMENTE PRESA
C COME SEED... DEVE ESSERE
C PRESA SOLO SE IL CLUSTER E' DOUBLE PEAKED...)
c------------------------------------------------------------------------
c 6 ottobre 2003
c Elena: CLSEEDCUT = 7 (old value 10)
c Elena: CLINCLCUT = 4 (old value 5)
iseed=-999 !cluster seed index initialization
do jl=1,nladders_view !1..3 !loops on ladders
first=1+nstrips_ladder*(jl-1) !1,1025,2049
last=nstrips_ladder*jl !1024,2048,3072
c X views have 1018 strips instead of 1024
if(mod(iv,2).eq.0) then
first=first+3
last=last-3
endif
do is=first,last !loop on strips in each ladder
if(is.le.iseed+1) goto 220
c-----------------------------------------
c after a cluster seed as been found,
c look for next one skipping one strip on the right
c (i.e. look for double peak cluster)
c-----------------------------------------
if(is.ne.first) then
if(value(is).le.value(is-1)) goto 220
endif
c-----------------------------------------
c skips cluster seed
c finding if strips values are descreasing (a strip
c can be a cluster seed only if previous strip value
c is lower)
c-----------------------------------------
if(value(is).gt.clseedcut(is)) then
ccc print*,"value(",is,")=",value(is),
ccc $ " .gt.clseedcut(",is,")=",clseedcut(is)
c-----------------------------------------
c possible SEED...
c-----------------------------------------
itemp=is
if(itemp.eq.last) goto 230 !estremo...
do while(value(itemp)
$ /sigma(iv,nvk(itemp),nst(itemp))
$ .le.value(itemp+1)
$ /sigma(iv,nvk(itemp+1),nst(itemp+1))) !BIAS: aggiustare il caso uguale!???
itemp=itemp+1
if(itemp.eq.last) goto 230 !stops if reaches last strip
enddo ! of the ladder
230 continue
c-----------------------------------------
c fownd SEED!!!
c-----------------------------------------
iseed=itemp
c----------------------------------------------------------
c after finding a cluster seed, checks also adjacent strips,
C and marks the ones exceeding clinclcut
c----------------------------------------------------------
ir=iseed !indici destro
il=iseed ! e sinistro
rmax=ir !estremo destro del cluster
lmax=il ! e sinistro
rstop=0 !initialize flags used to exit from
lstop=0 ! inclusion loop
do while(lstop.eq.0.or.rstop.eq.0) !shifts left and right from
ir=ir+1 !position index for strips on right side of
! cluster seed
il=il-1 !and for left side
c------------------------------------------------------------------------
c checks for last or first strip of the ladder
c------------------------------------------------------------------------
if(ir.gt.last) then !when index goes beyond last strip
rstop=1 ! of the ladder, change rstop flag in order
! to "help" exiting from loop
endif
if(il.lt.first) then !idem when index goes beyond
lstop=1 ! first strip of the ladder
endif
c------------------------------------------------------------------------
c check for clusters including more than nclstrp strips
c------------------------------------------------------------------------
if((rmax-lmax+1).ge.nclstrp) then
goto 210 !exits inclusion loop:
! lmax and rmax maintain last value
! NB .ge.!???
endif
c------------------------------------------------------------------------
c marks strips exceeding inclusion cut
c------------------------------------------------------------------------
if(rstop.eq.0) then !if last strip of the ladder or last
! over-cut strip has not been reached
if(value(ir).gt.clinclcut(ir)) then !puts in rmax the
rmax=ir ! last right over-cut strip
else
rstop=1 !otherwise cluster ends on right and rstop
endif ! flag=1 signals it
endif
if(lstop.eq.0) then
if(value(il).gt.clinclcut(il)) then
lmax=il
else
lstop=1
endif
endif
enddo !ends strip inclusion loop
210 continue !jumps here if more than nclstrp have been included
multtemp=rmax-lmax+1 !stores multiplicity in temp
! variable. NB rmax and lmax can change later in
! order to include enough strips to calculate eta3
! and eta4. so mult is not always equal to cllength
c------------------------------------------------------------------------
c NB per essere sicuro di poter calcolare eta3 e eta4 devo includere
c sempre e comunque le 2 strip adiacenti al cluster seed e quella
c adiacente ulteriore dalla parte della piu' alta fra queste due
c (vedi oltre...)!???
c------------------------------------------------------------------------
c nel caso di estremi del ladder...!???
c ho meno di 4 strip nel cluster --> se sono sui bordi o quasi del ladder
c costruisco il cluster ad hoc e poi esco, se non sono sui bordi o quasi
c vado oltre (aggiungero' quindi strip a sx e dx in modo da poter calcolare
c eta3e4)
if((rmax-lmax+1).lt.4) then
if(iseed.eq.first) then !estremi...
rmax=iseed+2 !NB in questo modo puo' anche capitare di
lmax=iseed ! includere strip sotto taglio di inclusione
goto 250 ! che non serviranno per eta3e4!???
endif
if(iseed.eq.last) then !estremi...
rmax=iseed
lmax=iseed-2 !NB 2 e non 3, perche' altrimenti sarei in
goto 250 ! ((rmax-lmax+1).lt.4).eq.false. !???
endif !NMB questo e' l'unico caso di cllength=3!???
if(iseed.eq.first+1) then !quasi estremi...
rmax=iseed+2
lmax=iseed-1
goto 250
endif
if(iseed.eq.last-1) then
rmax=iseed+1
lmax=iseed-2
goto 250
endif
c se ho 4 o piu' strip --> se sono sui bordi esco, se sono sui quasi bordi
c includo la strip del bordo
else
if(iseed.eq.first) goto 250 !estremi... non includo altro
if(iseed.eq.last) goto 250
if(iseed.eq.first+1) then !quasi estremi... mi assicuro di
lmax=first ! avere le strip adiacenti al seed
if((rmax-lmax+1).gt.nclstrp) rmax=rmax-1 !NB effetto
goto 250 ! coperta: se la lunghezza del cluster era gia'
endif ! al limite (nclstrp), per poter aggiungere questa
! strip a sinistra devo toglierne una a destra...!???
if(iseed.eq.last-1) then
rmax=last
if((rmax-lmax+1).gt.nclstrp) lmax=lmax+1
goto 250
endif
endif
c------------------------------------------------------------------------
c be sure to include in the cluster the cluster seed with its 2 adjacent
c strips, and the one adjacent to the greatest between this two strip, as the
c fourth one. if the strips have the same value (!) the fourth one is chosen
c as the one having the greatest value between the second neighbors
c------------------------------------------------------------------------
if(value(iseed+1).eq.value(iseed-1)) then
if(value(iseed+2).ge.value(iseed-2)) then !??? qui cmq c'e'
diff=(iseed+2)-rmax
if(diff.gt.0) then
rmax=rmax+diff
if((rmax-lmax+1).gt.nclstrp) then
lmax=rmax-nclstrp+1
endif
endif
diff=(iseed-1)-lmax
if(diff.lt.0) then
lmax=lmax+diff
if((rmax-lmax+1).gt.nclstrp) then
rmax=lmax+nclstrp-1
endif
endif
else
diff=(iseed-2)-lmax
if(diff.lt.0) then
lmax=lmax+diff
if((rmax-lmax+1).gt.nclstrp) then
rmax=lmax+nclstrp-1
endif
endif
diff=(iseed+1)-rmax
if(diff.gt.0) then
rmax=rmax+diff
if((rmax-lmax+1).gt.nclstrp) then
lmax=rmax-nclstrp+1
endif
endif
endif
elseif(value(iseed+1).gt.value(iseed-1)) then
c !??? sposto il limite del cluster a destra per includere sempre le strip
c necessarie al calcolo di eta-i
c se il cluster diventa troppo lungo lo accorcio a sinistra per avere non piu'
c di nclstrp (in questo caso sono sicuro di aver gia' incluso le strip
c necessarie al calcolo di eta-i a sinistra, quindi se voglio posso uscire)
diff=(iseed+2)-rmax
if(diff.gt.0) then
rmax=rmax+diff
if((rmax-lmax+1).gt.nclstrp) then
lmax=rmax-nclstrp+1
c goto 250
endif
endif
diff=(iseed-1)-lmax
if(diff.lt.0) then
lmax=lmax+diff
if((rmax-lmax+1).gt.nclstrp) then
rmax=lmax+nclstrp-1
c goto 250 !inutile!???
endif
endif
else
diff=(iseed-2)-lmax
if(diff.lt.0) then
lmax=lmax+diff
if((rmax-lmax+1).gt.nclstrp) then
rmax=lmax+nclstrp-1
c goto 250
endif
endif
diff=(iseed+1)-rmax
if(diff.gt.0) then
rmax=rmax+diff
if((rmax-lmax+1).gt.nclstrp) then
lmax=rmax-nclstrp+1
c goto 250 !inutile!???
endif
endif
endif
250 continue
c--------------------------------------------------------
c fills ntuple variables
c--------------------------------------------------------
nclstr1=nclstr1+1 !cluster number
ccc print*,nclstr1,multtemp
if(nclstr1.gt.nclstrmax) then !too many clusters for the event:
good1=0 ! event
nclstr1=0
totCLlength=0
flag_shower = .true.
if(verbose)print*,'Event ',eventn(1),
$ ': more than ',nclstrmax,' clusters'
goto 2000
endif
view(nclstr1)=iv !vista del cluster
ladder(nclstr1)=nld(iseed,iv) !ladder a cui appartiene il cluster seed
maxs(nclstr1)=iseed !strip del cluster seed
mult(nclstr1)=multtemp !molteplicita'
indstart(nclstr1)=ind !posizione dell'inizio del cluster nell'
! array clsignal
indmax(nclstr1)=indstart(nclstr1)+(iseed-lmax) !posizione del
! cluster seed nell'array clsignal
CLlength=rmax-lmax+1 !numero di strip del cluster
totCLlength=totCLlength+CLlength
dedx(nclstr1)=0
do j=lmax,rmax !stores sequentially cluter strip values in
clsignal(ind)=value(j) ! clsignal array
ind=ind+1
c if(value(j).gt.0)
if(value(j).gt.clinclcut(j))
$ dedx(nclstr1)=dedx(nclstr1)+value(j) !cluster charge
enddo
c--------------------------------------------------------
c
c--------------------------------------------------------
endif !end possible seed conditio
220 continue !jumps here to skip strips left of last seed
enddo ! end loop on strips
enddo !end loop on ladders
2000 continue
return
end
*---***---***---***---***---***---***---***---***
*
*
*
*
*
*---***---***---***---***---***---***---***---***
subroutine stripmask
* this routine set va1 and single-strip masks,
* on the basis of the VA1 mask saved in the DB
*
* mask(nviews,nva1_view,nstrips_va1) !strip mask
* mask_vk(nviews,nva1_view) !VA1 mask
*
include 'commontracker.f'
include 'level1.f'
include 'calib.f'
* init mask
do iv=1,nviews
do ivk=1,nva1_view
do is=1,nstrips_va1
mask(iv,ivk,is) = mask_vk(iv,ivk)
enddo
enddo
enddo
return
end