*************************************************************************
*
* Subroutine cncomp.f
*
* iterates common noise computation subroutine (./cnoise.f) and cluster
* cutting subroutine (./cutcn.f) till no more clusters are found
*
*************************************************************************
subroutine cncomp(i,j) !(view, VA1)
include 'commontracker.f'
include 'common_reduction.f'
include 'calib.f'
integer errflag !error flag to mark no signal free VA1
integer clstr_old(nstrips_va1) !flag storage vector
real signal(nstrips_va1) !"signal" (=adc-ped) value storage vector
real smean, ssigma !"signal" mean and sigma
real cut !"strange" strip exclusion cut
integer newclstr !flag to warn about new found clusters to be
! excluded from common noise computation
c call HBOOK1(20000+100*i+j,' ',30,0.,30.,0.) !???
c------------------------------------------------------------------------
c
c variables initialization
c
c------------------------------------------------------------------------
do k=1,nstrips_va1 !loops on strips
clstr(i,j,k)=1 !initializes signal affected strips flag
clstr_old(k)=1 !initializes signal affected strips storage
strange(i,j,k)=1 !initializes unusually high or low signal
enddo ! affected strips flag
newclstr=1 !flag to warn about new found signal
! affected strips
c------------------------------------------------------------------------
c
c high or low signal affected strips exclusion: computes "signal" (=adc-ped)
c mean value and sigma, and cuts from common noise computation strips
c whose ABS(signal) exceeds scut*sigma
c
c------------------------------------------------------------------------
countme=0 !???
666 continue !???
smean=0. !initialization
ssigma=0.
nstr=0
do k=1,nstrips_va1
nstr=nstr+strange(i,j,k) !uses only
if(mod(i,2).eq.1) then !odd strip ---> Y view
signal(k)= - (DBLE(adc(i,j,k))-pedestal(i,j,k)) !negative signal
else !even strip ---> X view
signal(k)= DBLE(adc(i,j,k))-pedestal(i,j,k) !positive signal
endif
smean=smean+signal(k)*strange(i,j,k)
ssigma=ssigma+(signal(k)**2)*strange(i,j,k)
c call HFILL(10000+100*i+j,signal(k),0.,1.) !???
enddo
smean=smean/nstr !strips value distribution mean
ssigma=SQRT((ssigma/nstr)-smean**2) !strips value distribution sigma
cut=scut*ssigma !exclusion cut
do k=1,nstrips_va1
c call HFILL(20000+100*i+j,ABS(signal(k)-smean)/ssigma,0.,1.) !???
if(ABS(signal(k)-smean).gt.cut) then
c print*,i,j,k,signal(k),abs(signal(k)),cut,strange(i,j,k) !???
strange(i,j,k)=0 !marks strips exceeding cut
endif
enddo ! in order not to use them in CN computation
countme=countme+1 !???
if (countme.le.3) goto 666 !???
c------------------------------------------------------------------------
c
c common noise computation
c
c------------------------------------------------------------------------
do while(newclstr.eq.1) !loops on this VA1 till no new signal
! affected strips are found
newclstr=0 !to exit from loop if no new cluster is
! found
errflag=0
call cnoise(i,j,errflag) !(view, VA1, error flag) computes common
! noise
c print*,cn(i,j) !???
if(errflag.eq.1) goto 10 !goes to next VA1: this one has no signal
! free strips...
call cutcn(i,j) !(view, VA1) excludes clusters from
! common noise calculation
ncs=0 !initializes number of strips not excluded by cncut
do k=1,nstrips_va1 !loops on strips
if(clstr(i,j,k).ne.clstr_old(k)) then !checks if there are
! new found clusters, and if so sets
newclstr=1 ! newclstr flag = 1
clstr_old(k)=clstr(i,j,k) !stores cluster flags in
endif ! clstr_old variable
iok=strange(i,j,k)*bad(i,j,k)*clstr(i,j,k)
ncs=ncs+iok !counts number of good strips for cn computation
enddo
enddo !ends do while loop when there are no new
! clusters
c call HFILL(666,FLOAT(ncs),0.,1.) !???
c$$$ if(ncs.lt.20) then !warns if too many strips have been excluded from CN
c$$$ ! computation
c$$$ print*,'cncomp: WARNING, LESS THAN 20 STRIPS PASSED CN CUT'
c$$$ $ //' ON VA1 ',j,', VIEW ',i !NB questo errore e' "un po'" in conflitto
c$$$ ! con quello che setta errflag (vedi cnoise.f)...
c$$$
c$$$ endif
10 continue
return
end
*************************************************************************
*
* Subroutine cnoise.f!DA COMMENTARE!???
*
* uses adc(nviews,nva1_view,nstrips_va1) and
* pedestal(nviews,nva1_view,nstrips_va1) variables to compute common noise,
* and fills cn(nviews,nva1_view) variable. in the computation only
* not-bad and not-signal-affected strips are used
* (bad(nviews,nva1_view,nstrips_va1) and
* clstr(nviews,nva1_view,nstrips_va1) flags)
*
* needs:
* - ./common_calib.f
*
* to be called inside ./cncomp.f
*
*************************************************************************
subroutine cnoise(i,j,gulp) !(view, VA1)
include 'commontracker.f'
include 'common_reduction.f'
include 'calib.f'
integer gulp !error flag
ncn=0 !number of strips in cn computation
cn(i,j)=0 !initializes cn variable
do k=1,nstrips_va1 !loops on strips
iok=strange(i,j,k)*bad(i,j,k)*clstr(i,j,k) !flag to mark strange, bad
! or signal affected strips
ccc print*,i,j,k,strange(i,j,k),bad(i,j,k),clstr(i,j,k),iok !???
cn(i,j)=cn(i,j) + (DBLE(adc(i,j,k)) - pedestal(i,j,k))*iok !sums ADC-PED
! values to compute common noise
ncn = ncn + iok !counts number of strips in cn computation
enddo
ccc print*,'ncn= ',ncn
if(ncn.eq.0) then !no signal free strips on this VA1...
print*,'cnoise: WARNING, NO SIGNAL FREE STRIPS ON VA1 ',j,
$ ', VIEW ',i
gulp=1
else
cn(i,j)=cn(i,j)/DBLE(ncn) !computes common noise
gulp=0 !resets error flag
endif
return
end
*************************************************************************
*
* Subroutine cutcn.f!DA COMMENTARE!???
*
* excludes strips with particle signals and/or noisy strips from common
* noise calculation, marking their clstr(nviews,nva1_view,nstrips_va1)
* flag:
* clstr=0 ---> not to be used in CN computation
* clstr=1 ---> to be used in CN computation
*
* needs:
* - ./common_calib.f
*
* to be called inside ./cncomp.f
*
*************************************************************************
subroutine cutcn(i,j) !(view, VA1)
include 'commontracker.f'
include 'common_reduction.f'
include 'calib.f'
integer skip !used to skip strips (see later...)
integer kr, kl !position indexes to check signal affected
! strips on right and left side of cluster
! seed
integer ir, il !flags to exit loop on reaching VA1 extremes
real valuec !cluster seed signal
real cut,stripcut !cluster seed cut
real valuel, valuer !left and right strips signal
real stripcnincut !strip include cut
skip = 0 !initializes skip
do k=1,nstrips_va1 !loops on strips searching for cluster seeds
if(k.le.skip) goto 20 !continues only if k strip has not been
! checked yet
clstr(i,j,k)=1 !reinitializes strip to be used in CN!???
! computation, in order to be able to exclude
! different strips at every CN computation loop
c------------------------------------------------------------------------
c
c selects cut according to view
c
c------------------------------------------------------------------------
if(mod(i,2).eq.1) then !odd strip ---> Y view
valuec= - (DBLE(adc(i,j,k))-cn(i,j)-pedestal(i,j,k)) !negative signal
cut=clcuty !sets Y cut to find cluster seeds
else !even strip ---> X view
valuec= DBLE(adc(i,j,k))-cn(i,j)-pedestal(i,j,k) !positive signal
cut=clcutx !sets X cut to find cluster seeds
endif
c------------------------------------------------------------------------
c
c seeks clusters
c
c------------------------------------------------------------------------
stripcut=cut*sigma(i,j,k) !cluster seed cut
c if(ABS(valuec).gt.stripcut) then !checks if signal exceeds threshold!???
if(valuec.gt.stripcut) then !checks if signal exceeds threshold
c$$$ print*,'cut',i,j,k,valuec,stripcut,adc(i,j,k),cn(i,j)
c$$$ $ ,pedestal(i,j,k) !???
clstr(i,j,k)=0 !if so, marks this strip as part of a
! cluster
c------------------------------------------------------------------------
c after finding a cluster seed, checks also adiacent strips, and marks
c the ones exceeding cnincut
c------------------------------------------------------------------------
kr=k !initializes position indexes to be equal to
kl=k ! cluster seed position
ir=0 !initialize flags used to exit from
il=0 ! inclusion loop
do while (il.eq.0.or.ir.eq.0) !shifts left and right from
! cluster seed till it finds a strip below
! the threshold, or till it reaches first or
! last VA1 strip
kr=kr+1 !position index for strips on right side of
! cluster seed
kl=kl-1 !and for left side
c------------------------------------------------------------------------
c checks for last or first strip
c------------------------------------------------------------------------
if(kr.gt.nstrips_va1.and.ir.eq.0) then !when index goes
ir=1 ! beyond last VA1 strip, change ir flag in
! order to "help" exiting from loop
skip=nstrips_va1+1 !sets skip beyond last strip: all
! strips on the right have been included in
! the cluster, so skips all next strips
! (goto 20 condition is now always true)
endif
if(kl.lt.1.and.il.eq.0) then !idem when index goes beyond
il=1 ! first strip
endif
c P.S.: the "....and.i#.eq.0" term in above conditions is needed. In
c fact, even if I reach a under-cut strip on the right (so I get ir=1),
c the "do while loop" continues till such strip will be found on the
c left too.
c Thus kl and kr (! too) keep increasing, and it can happen kr gets
c greater than nstrips_va1 before kl reaches a under-cut strip. In this
c case it would pass this "if condition", so setting skip=nstrips_va1+1
c and skipping right strips never checked, if the "....and.i#.eq.0" term
c weren't the: instead, including this part it won't pass it
c because when I found reach the last VA1 strip on the right I set ir=1.
c (AAAAAAHHHHHHHHH!!!!!!!!!!!)
c------------------------------------------------------------------------
c marks strips exceeding inclusion cut
c------------------------------------------------------------------------
c for right strips (kr index)
if(ir.eq.0) then !if last strip or last over-cut strip has
! not been reached
if(mod(i,2).eq.1) then !Y view
valuer= - (DBLE(adc(i,j,kr))-cn(i,j) !puts in valuer
$ -pedestal(i,j,kr)) ! right strip value
else !X view
valuer=DBLE(adc(i,j,kr))-cn(i,j)-pedestal(i,j,kr)
endif
stripcnincut=cnincut*sigma(i,j,kr) !defines include cut
c if(ABS(valuer).gt.stripcnincut) then !marks right strip if it !???
if(valuer.gt.stripcnincut) then !marks right strip if it
clstr(i,j,kr)=0 !exceedes include cut
c$$$ print*,'inclcut_r',i,j,kr,valuer,stripcnincut
c$$$ $ ,adc(i,j,kr),cn(i,j),pedestal(i,j,kr) !???
else
ir=1 !otherwise cluster ends and ir flag =1
! signals it
skip=kr !putting skip=kr, next k to be checked is
! k=kr
endif
endif
c for left strips (kl index)
if(il.eq.0) then !if first strip or last over-cut strip has
! not been reached
if (mod(i,2).eq.1) then !Y view
valuel= - (DBLE(adc(i,j,kl))-cn(i,j) !puts in valuel
$ -pedestal(i,j,kl)) ! left strip value
else !X view
valuel=DBLE(adc(i,j,kl))-cn(i,j)-pedestal(i,j,kl)
endif
stripcnincut=cnincut*sigma(i,j,kl) !defines include cut
c if(ABS(valuel).gt.stripcnincut) then !marks left strip if it!???
if(valuel.gt.stripcnincut) then !marks left strip if it
clstr(i,j,kl)=0 !exceedes include cut
c$$$ print*,'inclcut_l',i,j,kl,valuel,stripcnincut
c$$$ $ ,adc(i,j,kl),cn(i,j),pedestal(i,j,kl) !???
else
il=1 !otherwise cluster ends and il flag =1
! signals it
endif
endif
enddo !ends lateral strips loop
endif !ends cluster seed condition
20 continue !comes here if next strip on the right has
! already been included in a cluster
enddo !ends principal strip loop
return
end