************************************************************************* * * 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 call init_level1 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 GOOD1(DSPnumber(iv)) = 2 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 GOOD1(DSPnumber(iv))=3 goto 18 endif c ------------------------ c DSP-counter jump c ------------------------ if( $ eventn_old(iv).ne.0.and. !first event in this file $ eventn(iv).ne.1.and. !first event in run $ good_old(DSPnumber(iv)).ne.0.and. !previous event corrupted $ .true.)then if(eventn(iv).ne.(eventn_old(iv)+1))then GOOD1(DSPnumber(iv))=4 goto 18 endif endif c ------------------------ 18 continue endif enddo ngood = 0 do iv = 1,nviews eventn_old(iv) = eventn(iv) good_old(iv) = good1(iv) ngood = ngood + good1(iv) enddo c if(ngood.ne.0)print*,'* WARNING * 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-------------------------------------------------- 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-------------------------------------------------- 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 c NBNBNBNBNB mask per la striscia 1 !!!!!!!! if(mask(iv,ik,1).eq.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 enddo 100 format(' * WARNING * Event ',i7,' view',i3,': VK MASK ',24i1) if(ima.ne.0.and.debug)write(*,100)eventn(1),iv $ ,(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 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)) 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) else fshower(iv) = 1 c GOOD1(DSPnumber(iv)) = 11 !AHAHAHAHA ORRORE!! GOOD1(iv) = 11 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 c$$$ if(ngood.ne.0)print*,'* WARNING * Event ',eventn(1) c$$$ $ ,':LEVEL1 event status: ' c$$$ $ ,(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 dedx(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 morder 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 morder 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 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 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 mult_view(nclstr_view) = rmax-lmax+1 rmax_view(nclstr_view) = rmax lmax_view(nclstr_view) = lmax 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 dedx(nclstr1) = 0 do j=lmax_view(ic),rmax_view(ic) !stores sequentially cluter strip values in clsignal(ind) = value(j) ! clsignal array 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)) $ dedx(nclstr1) = dedx(nclstr1) + value(j) !cluster charge enddo c print*,'view ',iv,' -- save_cluster -- nclstr1: ' c $ ,nclstr1,maxs(nclstr1),mult(nclstr1),dedx(nclstr1) 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( 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 else mask(iv,ivk,is) = -1 $ * mask_vk(iv,ivk) !from DB $ * mask_vk_ev(iv,ivk) !from CN endif enddo return end