************************************************************************* * * 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,errflag) !(view, VA1) include 'commontracker.f' include 'level1.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 c ! excluded from common noise computation 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 c------------------------------------------------------------------------ c (september 2007) c remove from CN computation the first and the last 3 channels of c each X view, becouse they ar not connected to any strip c------------------------------------------------------------------------ if(mod(i,2).eq.0)then if(j.eq.1)then do k=1,3 strange(i,j,k)=0 enddo elseif(j.eq.nva1_ladder)then do k=nstrips_va1,nstrips_va1-2,-1 strange(i,j,k)=0 enddo endif endif newclstr=1 !flag to warn about new found signal c ! 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 ! ---> Y view signal(k) = - (DBLE(adc(i,j,k)) - pedestal(i,j,k)) !negative signal else ! ---> 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) enddo smean=smean/nstr !strips value distribution mean ssigma=SQRT((ssigma/nstr)-smean**2) !strips value distribution sigma cut=scut*ssigma !exclusion cut nco=0 nbo=0 do k=1,nstrips_va1 if(ABS(signal(k)-smean).gt.cut) then strange(i,j,k)=0 !marks strips exceeding cut c print*,i,j,k,signal(k),smean endif nco=nco+strange(i,j,k) nbo=nbo+bad(i,j,k) enddo ! in order not to use them in CN computation c$$$ if(i.eq.12.and.(j.eq.2.or.j.eq.3))then c$$$ print*,'view ',i,' vk ',j c$$$ print*,'ADC (1-51-128) = ',adc(i,j,1),adc(i,j,52),adc(i,j,128) c$$$ print*,' = ',smean c$$$ print*,'s = ',ssigma c$$$ print*,'nstrange = ',128-nco c$$$ print*,'nbad = ',128-nbo c$$$ endif countme = countme + 1 !??? if (countme.le.3) goto 666 !??? c------------------------------------------------------------------------ c c common noise computation c c----------------------------------------------------------------------- * loops on this VA1 till no new signal affected strips are found do while(newclstr.eq.1) newclstr=0 !to exit from loop if no new cluster is found errflag=0 call cnoise(i,j,errflag) !(view, VA1, error flag) computes cn 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 cn computation ncs=0 !initializes number of strips not excluded by cncut do k=1,nstrips_va1 !loops on strips * checks if there are new found clusters, and if so sets if(clstr(i,j,k).ne.clstr_old(k)) then newclstr=1 clstr_old(k)=clstr(i,j,k) endif 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 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 'level0.f' include 'level1.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 cnrms(i,j)=0 !initializes cn rms cnn(i,j)=0 !initialize cn flag do k=1,nstrips_va1 !loops on strips * tags strange, bad or signal-affected strips iok = strange(i,j,k)*bad(i,j,k)*clstr(i,j,k) cn(i,j) = cn(i,j) + (DBLE(adc(i,j,k)) - pedestal(i,j,k))*iok cnrms(i,j) = cnrms(i,j) $ + (DBLE(adc(i,j,k)) - pedestal(i,j,k)) $ *(DBLE(adc(i,j,k)) - pedestal(i,j,k))*iok ncn = ncn + iok !counts number of strips in cn computation enddo if(ncn.lt.NSTRIPMIN) then !no signal free strips on this VA1... if(ncn.eq.0)then if(debug.eq.1)print*,' WARNING - cnoise: ', $ 'no strips for CN computation on VA1 ',j, $ ', VIEW ',i,' >>> FAILED ' else if(debug.eq.1)print*,' WARNING - cnoise: ', $ 'less than ',NSTRIPMIN $ ,' strips for CN computation on VA1 ',j, $ ', VIEW ',i,' >>> FAILED ' endif gulp=1 cnn(i,j) = 0 else cn(i,j)=cn(i,j)/DBLE(ncn) !<<<< computes common noise cnrms(i,j)= SQRT( cnrms(i,j)/DBLE(ncn) - cn(i,j)**2 ) cnn(i,j) = ncn gulp=0 c$$$ print*,'Event ',eventn(1) c$$$ $ ,': cn(',i,',',j,')= ',cn(i,j),' ncn ',ncn if(debug.eq.1.and.ABS(cn(i,j)).gt.1000) $ print*,'Event ',eventn(1) $ ,': cn(',i,',',j,')= ',cn(i,j),' ncn ',ncn 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 'level1.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