************************************************************************* * * Program reduction.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 * * * * 11/9/2005 modified by david fedele to include general variables * ************************************************************************* program reduction include '../common/commontracker.f' include '../common/level0.f' include '../common/level1.f' include '../common/common_reduction.f' include '../common/calib.f' c------------------------------------------------------------------------ c c local variables c c------------------------------------------------------------------------ parameter(ncalibmax=50) c local variables character*74 data_file !data file name character*74 data_dir !data file name character*74 data_file_calib !calibration list file character*40 file_calib(ncalibmax) character*74 data_file_level0 character*74 data_file_level1 integer minevent !first event to be analysed integer nevent !number of events to be analysed character*24 processing_date integer which_calib_last parameter (lun_data_file=70) !data file id number parameter (lun_data_level0=71) !data file id number parameter (lun_data_level1=72) !data file id number parameter(lun_calib_list=66) !calibration list file id integer ind real value(nstrips_view) !per trovare i cluster real clseedcut(nstrips_view) real clinclcut(nstrips_view) common/searchcluster/value,clseedcut,clinclcut,ind integer nshowers logical flag_shower common/shower/nshowers,flag_shower COMMON/QUEST/IQUEST(100) !permette di ottenere ntuple funzionanti nonostante ! il messaggio dei 64K di RZOUT...!??? c***************************************************** cccccc 11/9/2005 modified by david fedele swcode1=202 c**************************************************** c------------------------------------------------------------------------ c c HBOOK initialization c c------------------------------------------------------------------------ call HLIMIT(NWPAWC) c------------------------------------------------------------------------ c c reads input informations (through < go_reduction) c c------------------------------------------------------------------------ call fdate(processing_date) write(*,101) $ processing_date 101 format(/ $ ,'*** *** *** *** *** *** *** *** *** *** *** *** ***',/ $ ,'* *',/ $ ,'* REDUCTION *',/ $ ,'* *',/ $ ,'*** *** *** *** *** *** *** *** *** *** *** *** ***',/ $ ,a24,/ $ ) 111 format(a) print*,'Data directory:' read(*,111)data_dir print*,data_dir print*,'File identifier: (DATE_NUM)' read(*,*)data_file print*,data_file c$$$ print*, c$$$ $ data_dir(1:LNBLNK(data_dir)) c$$$ $ //'DW_'//data_file(1:LNBLNK(data_file)) minevent=1 print*,'Maximum number of events to be analized:' !20000 read(*,*) nevent print*,nevent print*,'---------------------------------------------------' c------------------------------------------------------------- c------------------------------------------------------------- c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c reads calibration file c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c------------------------------------------------------------- c------------------------------------------------------------- print*,' ' print*,'OPENING CALIBRATION-LIST FILE:' 501 format(a,'DW_',a,'_calib.txt') write(data_file_calib,501) $ data_dir(1:LNBLNK(data_dir)) $ ,data_file(1:LNBLNK(data_file)) print*,data_file_calib open(lun_calib_list, $ FILE=data_file_calib(1:LNBLNK(data_file_calib)) $ ,STATUS='UNKNOWN' $ ,IOSTAT=iostat $ ) 113 format(i5,' ',a25) do i=1,ncalibmax read(lun_calib_list,113,IOSTAT=iostat) n_cal_list $ ,file_calib(i)(1:LNBLNK(file_calib(i))) if(iostat.ne.0)then ncal=i-1 goto 2000 endif print*,n_cal_list,' - ' $ ,file_calib(i)(1:LNBLNK(file_calib(i))) enddo 2000 close(lun_calib_list) which_calib_last=0 c-------------------------------------------------------------- c-------------------------------------------------------------- c c opening LEVEL 0 ntuple c c-------------------------------------------------------------- c-------------------------------------------------------------- 504 format(a,'DW_',a,'_level0.rz') write(data_file_level0,504) $ data_dir(1:LNBLNK(data_dir)) $ ,data_file(1:LNBLNK(data_file)) print*,'' print*,'OPENING LEVEL0 FILE:' print*,data_file_level0 IQUEST(10)=65000 !permette di ottenere ntuple funzionanti nonostante ! il messaggio dei 64K di RZOUT...!??? call HROPEN(lun_data_level0, $ 'LEVEL0',data_file_level0,'QP',4096,istat) !opens rz if(istat.ne.0) goto 19 call HRIN(0,9999,0) call access_level0 c call HPRNTU(ntp_level0) call HNOENT(ntp_level0,iemax0) print*,' events',iemax0 c------------------------------------------------------------------------ c c LEVEL 1 ntuple booking c c------------------------------------------------------------------------ 503 format(a,'DW_',a,'_level1.rz') write(data_file_level1,503) $ data_dir(1:LNBLNK(data_dir)) $ ,data_file(1:LNBLNK(data_file)) print*,'' print*,'------------------------------------' print*,' Creating LEVEL1 rz file' print*,' ',data_file_level1 print*,'------------------------------------' C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C largest RZ file: IQUEST(10) records x LREC words x 4 byte C with the following settings: 65000 x 4096 x 4 = 1G C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - IQUEST(10)=65000 c permette di ottenere ntuple funzionanti nonostante c il messaggio dei 64K di RZOUT call HROPEN(lun_data_level1, $ 'LEVEL1',data_file_level1,'QNP',4096,istat) !opens rz if(istat.ne.0) goto 19 call book_level1 * ------------------------------------------------------- * STRIP MASK * ------------------------------------------------------- call stripmask(data_file) nshowers=0 nocalib=0 c counts number of events with more than nclstrmax clusters maxevent=minevent+nevent-1 nevent = 0 nevent_good = 0 do iev = minevent,MIN(iemax0,maxevent) !loop on events nevent = nevent + 1 if(mod(nevent,100).eq.0)print*,nevent call init_level1 call HCDIR('//LEVEL0',' ') call HGNT(ntp_level0,iev,ierr) !read event at LEVEL0 if(ierr.ne.0)then print*,'HGNT ---> error reading event ',iev good1 = .false. nev1 = nev1 + 1 goto 200 !jump to fill n-tuple endif c if(mod(iev,100).eq.0) c $ print*,' event number: ',iev C--------------------------------------------------- c***************************************************** cccccc 11/9/2005 modified by david fedele ccC variables in blocks EVENT and CPU are anyway filled C variables in blocks GENERAL and CPU are anyway filled c******************************************************** C in order to mantain sincronization among C events at different levels C--------------------------------------------------- c***************************************************** cccccc 11/9/2005 modified by david fedele c good1=good0 c nev1=nev0 c pkt_type1=pkt_type c pkt_num1=pkt_num c obt1=obt c which_calib1=which_calib good1=good0 nev1=nev0 which_calib1=which_calib pkt_type1=pkt_type pkt_num1=pkt_num obt1=obt cpu_crc1=cpu_crc do iv=1,12 crc1(iv)=crc(iv) enddo c******************************************************* C--------------------------------------------------- C retrieve info about calib to be used C--------------------------------------------------- if(which_calib.ne.which_calib_last.and. $ which_calib.ne.0)then data_file_calib= $ data_dir(1:LNBLNK(data_dir))// $ file_calib(which_calib) $ (1:LNBLNK(file_calib(which_calib))) c print*,data_file_calib print*,'' print*, $ '@ event ',nev0 $ ,' (CPU pkt N.',pkt_num,')' print*,'--> ',data_file_calib IQUEST(10)=65000 call HROPEN(lun_data_file, $ 'CALIB',data_file_calib,'QP',4096,istat) !opens if(istat.ne.0) goto 19 call HRIN(0,9999,0) call fillpedsig do iview=1,nviews call HDELET(id_hi_bad+iview) call HDELET(id_hi_ped+iview) call HDELET(id_hi_sig+iview) enddo call HREND('CALIB') close(lun_data_file) which_calib_last=which_calib elseif(which_calib.eq.0)then nocalib=nocalib+1 good1=.false. endif if(.not.good1) goto 200 !jump to fill n-tuple 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=.false. print*,'event ',iev,' >>>>> decode ERROR' goto 200 endif ccc /////////////////////////////////////////////////// da qui 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 if(mask_vk(iv,ik).eq.1)call cncomp(iv,ik) enddo enddo 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)) 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=.false. nev1 = nev0 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 enddo enddo nevent_good = nevent_good + 1 200 continue C--------------------------------------------- C come here if GOOD1=.false. C or the event has too many clusters C--------------------------------------------- call HCDIR('//LEVEL1',' ') call HFNT(ntp_level1) !fills LEVEL1 ntuple ccc /////////////////////////////////////////////////// a qui enddo !end loop on events c$$$ call HPRNT(ntp_level1) write(*,108) $ nevent $ ,nevent_good $ ,nshowers $ ,nclstrmax $ ,nocalib 108 format(/ $ ,'|------------------------------------------------|',/ $ ,'| Number of processed events |',/ $ ,'| Total ',i8,' |'/ $ ,'| Good ',i8,' |'/ $ ,'| (',i5,' events with > ',i5,' clusters) |'/ $ ,'| (',i5,' events without calibration) |'/ $ ,'|------------------------------------------------|',/ $ ) goto 9000 c------------------------------------------------------------------------ c c no error exit c c------------------------------------------------------------------------ print*,' ' print*,'REDUCTION SUCCESSFULLY COMPLETED' print*,' ' print*,' ' goto 9000 !happy ending c------------------------------------------------------------------------ c c data file opening error c c------------------------------------------------------------------------ 19 continue print*,' ' print*,'reduction: ERROR OPENING DATA FILE: ',data_file print*,' ' print*,' ' goto 9000 !the end c------------------------------------------------------------------------ c c data file event error c c------------------------------------------------------------------------ 20 continue print*,' ' print*,'reduction: ERROR IN DATA FILE: ',data_file print*,'NO EVENTS FOUND' print*,' ' print*,' ' goto 9000 !the end c------------------------------------------------------------------------ c c data file eleboration level error c c------------------------------------------------------------------------ 21 continue print*,' ' print*,'reduction: ERROR IN DATA FILE: ',data_file print*,'LEVEL0= ',level0 print*,' ' print*,' ' goto 9000 !the end c------------------------------------------------------------------------ c c level0 ntuple event reading error c c------------------------------------------------------------------------ 22 continue print*,' ' print*,'reduction: ERROR WHILE READING LEVEL0 NTUPLE, AT EVENT: ', $ iev print*,ierr print*,' ' print*,' ' goto 9000 !the end c------------------------------------------------------------------------ c c closes files and exits c c------------------------------------------------------------------------ 9000 continue call HCDIR('//LEVEL1',' ') call RZPURG(-1) call HROUT(ntp_level1,ICYCLE,'T') call HREND('level1') close(lun_data_level1) call HCDIR('//LEVEL0',' ') call HREND('level0') close(lun_data_level0) c$$$ close(lun_data_runinfo) stop end ***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...*** * * * * * * * * * ***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...*** subroutine book_level1 include '../common/commontracker.f' include '../common/level1.f' character*35 block1,block2,block3, $ block4,block5,block6,block7 $ ,block8 * LEVEL1 ntuple: call HBNT(ntp_level1,'LEVEL1',' ') c***************************************************** cccccc 11/9/2005 modified by david fedele c call HBNAME(ntp_level1,'EVENT',good1,'GOOD1:L,NEV1:I') c call HBNAME(ntp_level1,'GENERAL',good1,'GOOD1:L,NEV1:I c $ ,WHIC_CALIB1:I,SWCODE1:I') cccccc 06/10/2005 modified by elena vannuccini call HBNAME(ntp_level1,'GENERAL',good1,'GOOD1:L,NEV1:I $ ,WHIC_CALIB1:I,SWCODE1:I,CRC1(12):L') c********************************************************* call HBNAME(ntp_level1,'CPU',pkt_type1 $ ,'PKT_TYPE1:I $ ,PKT_NUM1:I $ ,OBT1:I'// c******************************************************** cccccc 11/9/2005 modified by david fedele c $ ,WHICH_CALIB1:I') $ ',CPU_CRC1:L') c******************************************************** 411 format('NCLSTR1:I::[0,',I4,']') !cluster number 412 format(',VIEW(NCLSTR1):I::[0,',I2,']') !view the cluster belongs to 413 format(',LADDER(NCLSTR1):I::[0,',I1,']') !ladder the cluster belongs to 414 format(',MAXS(NCLSTR1):I::[0,',I5,']') !cluster stip carrying the largest ! signal value 415 format(',MULT(NCLSTR1):I::[0,',I2,']') !cluster multiplicity 416 format(',INDSTART(NCLSTR1):I::[0,',I5,']') !cluster starting point index in ! clsignal array 417 format(',INDMAX(NCLSTR1):I::[0,',I5,']') !cluster maximum point index in $ ! clsignal array 418 format(',TOTCLLENGTH:I::[0,',I5,']') !sum of all clusters length write(block1,411) nclstrmax !maximum number of clusters per event write(block2,412) nviews !number of views write(block3,413) nladders_view !number of ladders per view write(block4,414) nstrips_view !number of strips per view write(block5,415) nclstrp !maximum number of strips in a cluster write(block6,416) maxlength !maximum number of strip belonging to clusters write(block7,417) maxlength ! for the whole event write(block8,418) maxlength call HBNAME(ntp_level1,'CLUSTER',nclstr1,block1//block2//block3 $ //block4//block5//',DEDX(NCLSTR1):r'//block6//block7 c***************************************************** cccccc 11/9/2005 modified by david fedele c $ //block8//',CLSIGNAL(TOTCLLENGTH):R') c $ //block8//',CLSIGNAL(TOTCLLENGTH):R,CRC1(12):L') cccccc 06/10/2005 modified by elena vannuccini $ //block8//',CLSIGNAL(TOTCLLENGTH):R') c***************************************************** 419 format('CNEV(',I4,',',I4,')') !common noise of the event for a certain view write(block8,419) nviews,nva1_view ! and VA1 call HBNAME(ntp_level1,'CNOISE',cnev,block8) return end c............................................................. subroutine init_level1 include '../common/commontracker.f' include '../common/level1.f' include '../common/level0.f' 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***************************************************** cccccc 11/9/2005 modified by david fedele crc1(iv)=.false. c************************************************** do ik=1,nva1_view cnev(iv,ik)=0 enddo enddo return end *---***---***---***---***---***---***---***---*** * * * * * *---***---***---***---***---***---***---***---*** subroutine search_cluster(iv) include '../common/commontracker.f' include '../common/level1.f' include '../common/calib.f' integer ind real value(nstrips_view) !per trovare i cluster real clseedcut(nstrips_view) real clinclcut(nstrips_view) common/searchcluster/value,clseedcut,clinclcut,ind integer nshowers logical flag_shower common/shower/nshowers,flag_shower 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 c 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 c----------------------------------------- c possible SEED... c----------------------------------------- itemp=is if(itemp.eq.last) goto 230 !estremo... c$$$ do while(value(itemp) c$$$ $ /sigma(iv,nvk(itemp),nst(itemp)) c$$$ $ .le.value(itemp+1) c$$$ $ /sigma(iv,nvk(itemp+1),nst(itemp+1))) !BIAS: aggiustare il caso uguale!??? do while(value(itemp) $ .le.value(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 if(nclstr1.gt.nclstrmax) then !too many clusters for the event: nshowers=nshowers+1 ! skips variables filling and go to next good1=.false. ! event nclstr1=0 totCLlength=0 flag_shower = .true. print*,'Event ',nev1, $ ': 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(data_file) * this routine set va1 and single-strip masks, * on the basis of date and id downlink * * mask(nviews,nva1_view,nstrips_va1) !strip mask * mask_vk(nviews,nva1_view) !VA1 mask * include '../common/commontracker.f' include '../common/level1.f' include '../common/calib.f' character*10 data_file character*3 aid character*6 adate integer id integer date * ---------------------- * retrieve date and id aid=data_file(8:10) adate=data_file(1:6) READ (aid, '(I3)'), id READ (adate, '(I6)'), date * ---------------------- * init mask do iv=1,nviews do ivk=1,nva1_view mask_vk(iv,ivk)=1 do is=1,nstrips_va1 mask(iv,ivk,is)=1 enddo enddo enddo * --------------------- * VIEW 2 - VK 23-24 * couple of vk damaged during integration if(date.ge.50208)then print*,'MASK: view 2 - vk 23/24' mask_vk(2,23)=0 mask_vk(2,24)=0 do is=1,nstrips_va1 mask(2,23,is)=0 mask(2,24,is)=0 enddo endif * --------------------- * VIEW 7 - VK 11-12 if(date.ge.50209)then if(.not.(date.eq.50209.and.id.le.6)) then print*,'MASK: view 7 - vk 11/12' mask_vk(7,11)=0 mask_vk(7,12)=0 do is=1,nstrips_va1 mask(7,11,is)=0 mask(7,12,is)=0 enddo endif endif * --------------------- * VIEW 7 - VK 21-22 if(date.ge.50316)then print*,'MASK: view 7 - vk 21/22' mask_vk(7,21)=0 mask_vk(7,22)=0 do is=1,nstrips_va1 mask(7,21,is)=0 mask(7,22,is)=0 enddo endif * --------------------- * VIEW 12 - VK 1-2-3-4 if((date.eq.50317).and.(id.le.3))then print*,'MASK: view 12 - vk 1/2/3/4' mask_vk(12,1)=0 mask_vk(12,2)=0 mask_vk(12,3)=0 mask_vk(12,4)=0 do is=1,nstrips_va1 mask(12,1,is)=0 mask(12,2,is)=0 mask(12,3,is)=0 mask(12,4,is)=0 enddo endif * --------------------- * VIEW 7 - VK 5-6 if(date.ge.50320)then if(.not.(date.eq.50320.and.id.le.3)) then print*,'MASK: view 7 - vk 5/6' mask_vk(7,5)=0 mask_vk(7,6)=0 do is=1,nstrips_va1 mask(7,5,is)=0 mask(7,6,is)=0 enddo endif endif * --------------------- * VIEW 7 - VK 13-14 if(date.ge.50320)then if(.not.(date.eq.50320.and.id.le.5)) then print*,'MASK: view 7 - vk 13/14' mask_vk(7,13)=0 mask_vk(7,14)=0 do is=1,nstrips_va1 mask(7,13,is)=0 mask(7,14,is)=0 enddo endif endif *** SAMARA *** SAMARA *** SAMARA * it needs further checks... * --------------------- * VIEW 7 - VK 9-10 * VIEW 12 - VK 1-2-3-4 if((date.eq.50516).and.(id.le.8))then print*,'MASK: view 7 - vk 9/10' print*,'MASK: view 12 - vk 1/2/3/4' mask_vk(7,9)=0 mask_vk(7,10)=0 mask_vk(12,1)=0 mask_vk(12,2)=0 mask_vk(12,3)=0 mask_vk(12,4)=0 do is=1,nstrips_va1 mask(7,9,is)=0 mask(7,10,is)=0 mask(12,1,is)=0 mask(12,2,is)=0 mask(12,3,is)=0 mask(12,4,is)=0 enddo endif * --------------------- * VIEW 7 - VK 9-10 if(date.ge.50516)then if(.not.(date.eq.50516.and.id.le.8)) then print*,'MASK: view 7 - vk 9/10' mask_vk(7,9)=0 mask_vk(7,10)=0 do is=1,nstrips_va1 mask(7,9,is)=0 mask(7,10,is)=0 enddo endif endif * --------------------- * VIEW 12 - VK 7-8 if(date.ge.50523)then if(.not.(date.eq.50523.and.id.le.3)) then print*,'MASK: view 12 - vk 7/8' mask_vk(12,7)=0 mask_vk(12,8)=0 do is=1,nstrips_va1 mask(12,7,is)=0 mask(12,8,is)=0 enddo endif endif return end