************************************************************************* * * 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' data eventn_old/nviews*0/ integer ierror ierror = 0 c$$$ debug = .true. c$$$ verbose = .true. c$$$ warning = .true. c$$$ print*,debug,verbose,warning c$$$ debug=1 c$$$ verbose=1 c$$$ warning=1 * ////////////////////////// * initialize some parameters * ////////////////////////// call init_level1 c debug=.true. if(debug.eq.1)print*,'-- check LEVEL0 status' ievco=-1 mismatch=0 c good1 = good0 c-------------------------------------------------- c check the LEVEL0 event status for missing c sections or DSP alarms c ==> set the variable GOOD1(12) c-------------------------------------------------- do iv=1,nviews if(DSPnumber(iv).gt.0.and.DSPnumber(iv).le.12)then c ------------------------ c GOOD c ------------------------ GOOD1(DSPnumber(iv))=0 !OK c ------------------------ c CRC error c ------------------------ if(crc(iv).eq.1) then c GOOD1(DSPnumber(iv)) = 2 c GOOD1(DSPnumber(iv)) = GOOD1(DSPnumber(iv)) + 2**1 GOOD1(DSPnumber(iv)) = ior(GOOD1(DSPnumber(iv)),2**1) 102 format(' * WARNING * Event ',i7,' view',i3 $ ,' CRC error') if(debug.eq.1)write(*,102)eventn(1),DSPnumber(iv) c goto 18 !next view endif c ------------------------ c online-software alarm c ------------------------ if( $ fl1(iv).ne.0.or. $ fl2(iv).ne.0.or. $ fl3(iv).ne.0.or. $ fl4(iv).ne.0.or. $ fl5(iv).ne.0.or. $ fl6(iv).ne.0.or. $ fc(iv).ne.0.or. $ DATAlength(iv).eq.0.or. $ .false.)then c GOOD1(DSPnumber(iv))=3 c GOOD1(DSPnumber(iv)) = GOOD1(DSPnumber(iv)) + 2**2 GOOD1(DSPnumber(iv)) = ior(GOOD1(DSPnumber(iv)),2**2) 103 format(' * WARNING * Event ',i7,' view',i3 $ ,' software alarm') if(debug.eq.1)write(*,103)eventn(1),DSPnumber(iv) c goto 18 endif c ------------------------ c DSP-counter jump c ------------------------ c commentato perche` non e` un controllo significativo nel caso in cui c la subroutine venga chiamata per riprocessare l'evento c sostituito con un check dei contatori dei vari dsp c$$$ if( c$$$ $ eventn_old(iv).ne.0.and. !first event in this file c$$$ $ eventn(iv).ne.1.and. !first event in run c$$$ $ good_old(DSPnumber(iv)).ne.0.and. !previous event corrupted c$$$ $ .true.)then c$$$ c$$$ if(eventn(iv).ne.(eventn_old(iv)+1))then c$$$c GOOD1(DSPnumber(iv))=4 c$$$c GOOD1(DSPnumber(iv)) = GOOD1(DSPnumber(iv)) + 2**3 c$$$ GOOD1(DSPnumber(iv)) = ior(GOOD1(DSPnumber(iv)),2**3) c$$$ 104 format(' * WARNING * Event ',i7,' view',i3 c$$$ $ ,' counter jump ',i10,i10) c$$$ if(debug.eq.1)write(*,104)eventn(1),DSPnumber(iv) c$$$ $ ,eventn_old(iv),eventn(iv)) c$$$ goto 18 c$$$ endif c$$$ c$$$ endif c ------------------------ c 18 continue c ------------------------ c DSP-counter c ------------------------ if( DSPnumber(iv).ne.0.and.GOOD1(DSPnumber(iv)).ne.1)then if(iv.ne.1.and.ievco.ne.-1)then if( eventn(iv).ne.ievco )then mismatch=1 endif endif ievco = eventn(iv) endif endif enddo c print*,'*** ',(eventn(iv),iv=1,12) if(mismatch.eq.1.and.debug.eq.1) $ print*,' * WARNING * DSP counter mismatch: ' $ ,(eventn(iv),iv=1,12) ngood = 0 do iv = 1,nviews if(mismatch.eq.1.and.GOOD1(iv).ne.1) $ GOOD1(iv)=ior(GOOD1(iv),2**3) eventn_old(iv) = eventn(iv) good_old(iv) = good1(iv) ngood = ngood + good1(iv) enddo c$$$ if(verbose.and.ngood.ne.0)print*,'* WARNING * Event ',eventn(1) c$$$ $ ,':LEVEL0 event status: ' c$$$ $ ,(good1(i),i=1,nviews) c-------------------------------------------------- c read the variable DATATRACKER from LEVEL0 c and fill the variable ADC (invertin view 11) c-------------------------------------------------- if(debug.eq.1)print*,'-- fill ADC vectors' call filladc(iflag) if(iflag.ne.0)then ierror = 220 endif c-------------------------------------------------- c computes common noise for each VA1 c (excluding strips with signal, c tagged with the flag CLSTR) c-------------------------------------------------- if(debug.eq.1)print*,'-- compute CN' do iv=1,nviews ima=0 do ik=1,nva1_view cn(iv,ik) = 0 cnrms(iv,ik) = 0 cnn(iv,ik) = -1 iflag = 0 mask_vk_ev(iv,ik) = 1 call stripmask(iv,ik) !compute mask(i,j,k), combining VA1-masks * -------------------------------------- * if chip is not masked ---> evaluate CN * -------------------------------------- if( mask(iv,ik,1).eq.1 ) then !!!NBNB mask per la striscia 1 !!!!!!!! call cncomp(iv,ik,iflag) if(iflag.ne.0)then ima=ima+1 mask_vk_ev(iv,ik)=0 ierror = 220 endif call stripmask(iv,ik) !compute mask(i,j,k), combining VA1-masks endif enddo 100 format(' * WARNING * Event ',i7,' view',i3,': VK MASK ',24i1) if(ima.ne.0.and.verbose.eq.1)write(*,100)eventn(1),iv $ ,(mask_vk_ev(iv,ik),ik=1,nva1_view) c if(ima.ne.0)write(*,100)eventn(1),iv c $ ,(mask_vk_ev(iv,ik),ik=1,nva1_view) enddo cc call stripmask !compute mask(i,j,k), combining mask_vk_ev and mask_vk c--------------------------------------------- c loops on views, VA1 and strips, c and computes strips signals using c badstrip, pedestals, and c sigma informations from histograms c--------------------------------------------- ind=1 !clsignal array index if(debug.eq.1)print*,'-- search clusters' 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 c print*,iv,nvk(is),nst(is),adc(iv,nvk(is),nst(is)) c $ ,cn(iv,nvk(is)) c $ ,pedestal(iv,nvk(is),nst(is)) 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)) sat(is)=0 if( adc(iv,nvk(is),nst(is)).lt.adc_saty )sat(is)=1 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)) sat(is)=0 if( adc(iv,nvk(is),nst(is)).gt.adc_satx )sat(is)=1 endif enddo !end loop on strips (1) call search_cluster(iv) if(.not.flag_shower)then call save_cluster(iv) if(debug.eq.1)print*,'view ',iv,' #clusters ', nclstr_view else fshower(iv) = 1 c GOOD1(DSPnumber(iv)) = 11 !AHAHAHAHA ORRORE!! c GOOD1(iv) = 11 c GOOD1(iv) = GOOD1(iv) + 2**5 GOOD1(iv) = ior(GOOD1(iv),2**5) 101 format(' * WARNING * Event ',i7,' view',i3 $ ,' #clusters > ',i5,' --> MASKED') if(verbose.eq.1)write(*,101)eventn(1),iv,nclstrmax_view 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 cnrmsev(iv,ik) = cnrms(iv,ik) !assigns computed CN to ntuple variables cnnev(iv,ik) = cnn(iv,ik) !assigns computed CN to ntuple variables enddo enddo C--------------------------------------------- C come here if GOOD1=0 C or the event has too many clusters C--------------------------------------------- 200 continue ngood = 0 do iv = 1,nviews ngood = ngood + good1(iv) enddo if(verbose.eq.1.and.ngood.ne.0) $ print*,'* WARNING * Event ',eventn(1) $ ,':LEVEL1 event status: ' $ ,(good1(i),i=1,nviews) c------------------------------------------------------------------------ c c closes files and exits c c------------------------------------------------------------------------ RETURN END ***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...*** * * * * * * * * * ***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...*** subroutine init_level1 include 'commontracker.f' include 'level1.f' include 'level0.f' c good1 = 0 do iv=1,12 good1(iv) = 1 !missing packet enddo 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 sgnl(ic) = 0 whichtrack(ic) = 0 !assigned @ level2 enddo do id=1,maxlength !??? clsignal(id) = 0. clsigma(id) = 0. cladc(id) = 0. clbad(id) = 0. enddo do iv=1,nviews c crc1(iv)=0 do ik=1,nva1_view cnev(iv,ik) = 0 cnnev(iv,ik) = 0 enddo fshower(iv) = 0 enddo return end *---***---***---***---***---***---***---***---*** * * * * * *---***---***---***---***---***---***---***---*** subroutine search_cluster(iv) include 'commontracker.f' include 'level0.f' include 'level1.f' include 'calib.f' include 'common_reduction.f' c local variables integer rmax,lmax !estremi del cluster integer rstop,lstop integer first,last integer fsat,lsat external nst iseed=-999 !cluster seed index initialization inext=-999 !index where to start new cluster search flag_shower = .false. nclstr_view=0 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 * 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 c--------------------------------------------- c new-cluster search starts at index inext c--------------------------------------------- if(is.lt.inext) goto 220 ! next strip if(value(is).gt.clseedcut(is)) then c----------------------------------------- c possible SEED... c----------------------------------------- itemp = is fsat = 0 ! first saturated strip lsat = 0 ! last saturated strip if(itemp.eq.last) goto 230 !estremo... c ------------------------ c search for first maximum c ------------------------ do while( $ value(itemp).le.value(itemp+1) $ .and.value(itemp+1).gt.clseedcut(itemp+1)) itemp = itemp+1 if(itemp.eq.last) goto 230 !stops if reaches last strip if(sat(itemp).eq.1) goto 230 !stop if reaches a saturated strip enddo ! of the ladder 230 continue c ----------------------------- c check if strips are saturated c ----------------------------- if( sat(itemp).eq.1 )then fsat = itemp lsat = itemp if(itemp.eq.last) goto 231 !estremo... do while( sat(itemp+1).eq.1 ) itemp = itemp+1 lsat = itemp if(itemp.eq.last) goto 231 !stops if reaches last strip enddo endif 231 continue c--------------------------------------------------------------------------- c fownd SEED!!! c (if there are saturated strips, the cluster is centered in the middle) c--------------------------------------------------------------------------- if(fsat.eq.0.and.lsat.eq.0)then iseed = itemp ! <<< SEED else iseed = int((lsat+fsat)/2) ! <<< SEED c$$$ print*,'saturated strips ',fsat,lsat,iseed c$$$ print*,'--> ',(value(ii),ii=fsat,lsat) endif c--------------------------------------------------------------- c after finding a cluster seed, checks also adjacent strips, C and tags 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 !index for right side il=il-1 !index for left side c------------------------------------------------------------------------ c checks for last or first strip of the ladder c------------------------------------------------------------------------ if( ir.gt.last ) rstop = 1 if( il.lt.first ) lstop = 1 c------------------------------------------------------------------------ c add strips exceeding inclusion cut c------------------------------------------------------------------------ if( (rmax-lmax+1).ge.nclstrp )goto 210 !exits inclusion loop if(rstop.eq.0) then !if right cluster border has not been reached if(value(ir).gt.clinclcut(ir)) then rmax=ir !include a strip on the right else rstop=1 !cluster right end endif endif if( (rmax-lmax+1).ge.nclstrp )goto 210 !exits inclusion loop if(lstop.eq.0) then !if left cluster border has not been reached if(value(il).gt.clinclcut(il)) then lmax=il !include a strip on the left else lstop=1 !cluster left end endif endif c if( (rmax-lmax+1).ge.nclstrp )goto 210 !exits inclusion loop enddo !ends strip inclusion loop goto 211 210 continue !jumps here if more than nclstrp have been included c print*,'>>> nclstrp! ' 211 continue c----------------------------------------- c end of inclusion loop! c----------------------------------------- c------------------------------------------------------------------------ c adjust the cluster in order to have at least a strip around the seed c------------------------------------------------------------------------ if(iseed.eq.lmax.and.lmax.ne.first)then lmax = lmax-1 if( (rmax-lmax+1).gt.nclstrp )rmax=rmax-1 endif if(iseed.eq.rmax.and.rmax.ne.last )then rmax = rmax+1 if( (rmax-lmax+1).gt.nclstrp )lmax=lmax+1 endif c------------------------------------------------------------------------------- c adjust the cluster in order to have at least ANOTHER strip around the seed c------------------------------------------------------------------------------- c$$$ if(iseed-1.eq.lmax.and.lmax.ne.first)then c$$$ lmax = lmax-1 c$$$ if( (rmax-lmax+1).gt.nclstrp )rmax=rmax-1 c$$$ endif c$$$ if(iseed+1.eq.rmax.and.rmax.ne.last )then c$$$ rmax = rmax+1 c$$$ if( (rmax-lmax+1).gt.nclstrp )lmax=lmax+1 c$$$ endif c--------------------------------------------------- c now we have 5 stored-strips around the maximum c--------------------------------------------------- c------------------------------------------------------------------------ c adjust the cluster in order to store a minimum number of strips c------------------------------------------------------------------------ do while( (rmax-lmax+1).lt.nclstrpmin ) vl = -99999 vr = -99999 if(lmax-1.ge.first) vl = value(lmax-1) if(rmax+1.le.last ) vr = value(rmax+1) if(vr.ge.vl) then rmax = rmax+1 else lmax = lmax-1 endif enddo c-------------------------------------------------------- c store cluster info c-------------------------------------------------------- nclstr_view = nclstr_view + 1 !cluster number if(nclstr_view.gt.nclstrmax_view) then !too many clusters for the view: c$$$ if(verbose) print*,'Event ',eventn(1), c$$$ $ ': more than ',nclstrmax_view c$$$ $ ,' clusters on view ',iv flag_shower = .true. goto 2000 endif ladder_view(nclstr_view) = nld(iseed,iv) maxs_view(nclstr_view) = iseed rmax_view(nclstr_view) = rmax lmax_view(nclstr_view) = lmax c mult_view(nclstr_view) = rmax-lmax+1 mult_view(nclstr_view) = 0 do ii=lmax,rmax if(value(ii).gt.clinclcut(ii)) $ mult_view(nclstr_view) = mult_view(nclstr_view)+1 enddo c$$$ if(rmax-lmax+1.gt.25) c$$$ $ print*,'view ',iv c$$$ $ ,' cl ',nclstr_view,' mult ',rmax-lmax+1 c------------------------------------------------------------------------ c search for a double peak inside the cluster c------------------------------------------------------------------------ inext = rmax+1 !<< index where to start new-cluster search vmax = 0 vmin = value(iseed) imax = iseed imin = iseed do iss = max(iseed+1,lsat+1),rmax if( value(iss).lt.vmin )then if( imax.ne.iseed )goto 221 !found dowble peek imin = iss vmin = value(iss) else delta = value(iss) - value(imin) cut = sqrt(clinclcut(iss)**2 + clinclcut(imin)**2) if( $ delta.gt.cut .and. $ value(iss).gt.clseedcut(iss).and. $ .true.)then if( value(iss).gt.vmax )then imax = iss vmax = value(iss) else goto 221 !found dowble peek endif endif endif enddo 221 continue if(imax.gt.iseed)then inext = imax !<< index where to start new-cluster search c$$$ print*,'--- double peek ---' c$$$ print*,(value(ii),ii=lmax,rmax) c$$$ print*,'seed ',iseed,' imin ',imin,' imax ',imax endif 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 save_cluster(iv) * * (080/2006 Elena Vannuccini) * Save the clusters view by view include 'commontracker.f' include 'level1.f' include 'calib.f' include 'common_reduction.f' integer CLlength !lunghezza in strip del cluster do ic=1,nclstr_view nclstr1 = nclstr1+1 view(nclstr1) = iv ladder(nclstr1) = ladder_view(ic) maxs(nclstr1) = maxs_view(ic) mult(nclstr1) = mult_view(ic) c posizione dell'inizio del cluster nell' array clsignal indstart(nclstr1) = ind c posizione del cluster seed nell'array clsignal indmax(nclstr1) = indstart(nclstr1) $ +( maxs_view(ic) - lmax_view(ic) ) CLlength = rmax_view(ic) - lmax_view(ic) + 1 !numero di strip salvate totCLlength = totCLlength + CLlength sgnl(nclstr1) = 0 do j=lmax_view(ic),rmax_view(ic) !stores sequentially cluter strip values in clsignal(ind) = value(j) ! clsignal array c$$$ print*,ind,clsignal(ind) ivk=nvk(j) ist=nst(j) clsigma(ind) = sigma(iv,ivk,ist) cladc(ind) = adc(iv,ivk,ist) clbad(ind) = bad(iv,ivk,ist) c clped(ind) = pedestal(iv,ivk,ist) ind=ind+1 c if(value(j).gt.0) if(value(j).gt.clinclcut(j)) $ sgnl(nclstr1) = sgnl(nclstr1) + value(j) !cluster charge enddo c$$$ print*,'view ',iv,' -- save_cluster -- nclstr1: ' c$$$ $ ,nclstr1,maxs(nclstr1),mult(nclstr1),sgnl(nclstr1) c$$$ print*,'----------------------' enddo return end *---***---***---***---***---***---***---***---*** * * * * * *---***---***---***---***---***---***---***---*** c$$$ subroutine stripmask c$$$ c$$$* this routine set va1 and single-strip masks, c$$$* on the basis of the VA1 mask saved in the DB c$$$* c$$$* mask(nviews,nva1_view,nstrips_va1) !strip mask c$$$* mask_vk(nviews,nva1_view) !VA1 mask c$$$* c$$$ include 'commontracker.f' c$$$ include 'level1.f' c$$$ include 'common_reduction.f' c$$$ include 'calib.f' c$$$ c$$$* init mask c$$$ do iv=1,nviews c$$$ do ivk=1,nva1_view c$$$ do is=1,nstrips_va1 c$$$c mask(iv,ivk,is) = mask_vk(iv,ivk) c$$$ if( mask_vk(iv,ivk) .ne. -1)then c$$$ mask(iv,ivk,is) = 1 c$$$ $ * mask_vk(iv,ivk) !from DB c$$$ $ * mask_vk_ev(iv,ivk) !from c$$$ $ * mask_vk_run(iv,ivk) !from CN c$$$ else c$$$ mask(iv,ivk,is) = -1 c$$$ $ * mask_vk(iv,ivk) !from DB c$$$ $ * mask_vk_ev(iv,ivk) !from CN c$$$ endif c$$$ enddo c$$$ enddo c$$$ enddo c$$$ c$$$ c$$$ return c$$$ end subroutine stripmask(iv,ivk) * ----------------------------------------------- * 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 'common_reduction.f' include 'calib.f' * init mask do is=1,nstrips_va1 * -------------------------------------------------------- * if VA1-mask from DB is 0 or 1, three masks are combined: * - from DB (a-priori mask) * - run-based (chip declared bad on the basis of ) * - event-based (failure in CN computation) * -------------------------------------------------------- c print*,iv,ivk c $ ,mask_vk(iv,ivk),mask_vk_ev(iv,ivk),mask_vk_run(iv,ivk) if( mask_vk(iv,ivk) .ne. -1)then mask(iv,ivk,is) = 1 $ * mask_vk(iv,ivk) !from DB $ * mask_vk_ev(iv,ivk) !from $ * mask_vk_run(iv,ivk) !from CN * ----------------------------------------------------------- * if VA1-mask from DB is -1 only event-based mask is applied * ----------------------------------------------------------- else mask(iv,ivk,is) = -1 $ * mask_vk(iv,ivk) !from DB $ * mask_vk_ev(iv,ivk) !from CN endif enddo return end