************************************************************************* * Program analysis.f * * - reads cluster information (LEVEL1, reduction.f output ntuple) * - perform track identification and fit * - create LEVEL2 ntuple * ************************************************************************* program momanhough include '../common/commontracker.f' include '../common/common_momanhough.f' c include '../common/common_level2debug.f' include '../common/common_mech.f' include '../common/common_xyzPAM.f' include '../common/common_mini_2.f' include '../common/calib.f' include '../common/level1.f' include '../common/level2.f' include '../common/momanhough_init.f' * flag to set debug mode logical DEBUG common/dbg/DEBUG logical DBG_FILLED data DBG_FILLED/.false./ * flag to chose PFA character*10 PFA common/FINALPFA/PFA c parameter (inf=1.e8) !just a huge number... * external functions external npl external acoordsi,coordsi,nld,coord,dcoord c------------------------------------------------------------------------ c c local variables c c------------------------------------------------------------------------ character*24 processing_date parameter (lun_data_level1=71) !data file id number parameter (lun_data_level2=72) !data file id number parameter (lun_data_calib=74) !data file id number character*74 data_file !data file name character*74 data_dir !data file name character*74 data_file_calib character*74 data_file_level1 character*74 data_file_level2 # ifndef TEST2003 parameter(ncalibmax=50) character*40 file_calib(ncalibmax) parameter(lun_calib_list=66) !calibration list file id integer which_calib_last # endif integer minevent !first event to be analysed logical FIMAGE ! COMMON/QUEST/IQUEST(100) c !permette di ottenere ntuple funzionanti nonostante c ! il messaggio dei 64K di RZOUT...!??? c***************************************************** cccccc 11/9/2005 modified by david fedele c swcode=202 cccccc 12/10/2005 modified by Elena Vannuccini swcode=300 c**************************************************** c------------------------------------------------------------------------ c c HBOOK initialization c c------------------------------------------------------------------------ call HLIMIT(NWPAWC) c------------------------------------------------------------------------ c c reads input informations c c------------------------------------------------------------------------ call fdate(processing_date) write(*,101) $ processing_date 101 format(/ $ ,'*** *** *** *** *** *** *** *** *** *** *** *** ***',/ $ ,'* *',/ $ ,'* ANALYSIS *',/ $ ,'* *',/ $ ,'*** *** *** *** *** *** *** *** *** *** *** *** ***',/ $ ,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 minevent=1 print*,'Maximum number of events to be analized:' !20000 read(*,*) ntotev print*,ntotev print*,'Position-finding algorythm: (GOG2,ETA2)' read(*,*)PFA PFA=PFA(1:lnblnk(PFA)) print*,PFA print*,'DEBUG mode: (T, F)' 11 format(l1) read(*,11)DEBUG print*,DEBUG print*,'---------------------------------------------------' ****** INITIALIZATIONS ************************************* * 1) read charge-correlation parameters print*,' ' print*,'- read charge-correlation parameters' print*,' ' call readchargeparam * 1) read mip parameters print*,' ' print*,'- read mip parameters' print*,' ' call readmipparam * 2) read z coordinates of the planes print*,' ' print*,'- read z coordinates of the planes' print*,' ' call mech_sensor !reads sensors centres coordinates do ip=1,nplanes fitz(ip)=z_mech_sensor(ip,1,1)*0.1 !cm * gets planes mechanical z positions * (in mm) and sets them in micrometers enddo * 3) read eta PFA parameters print*,' ' print*,'- read P.F.A. parameters' print*,' ' call readetaparam print*,' ' print*,'- First guess P.F.A. >>>> ',PFAdef print*,' ' print*,'- Final P.F.A. >>>> ',PFA print*,' ' * 4) read magnetic field map print*,' ' print*,'- read magnetic field map' print*,' ' c call read_B(5) !legge il nome da STI!!.. temporaneo c call read_B_2maps(5) !legge il nome da STI!!.. temporaneo call read_B * 5) read allignment parameters print*,' ' print*,'- read aligment parameters' print*,' ' call readalignparam ************************************************************ c------------------------------------------------------------------------ c c LEVEL 2 ntuple booking c c------------------------------------------------------------------------ 503 format(a,'DW_',a,'_level2.rz') write(data_file_level2,503) $ data_dir(1:LNBLNK(data_dir)) $ ,data_file(1:LNBLNK(data_file)) print*,'' print*,'------------------------------------' print*,' Creating LEVEL2 rz file' print*,' ',data_file_level2 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_level2, $ 'LEVEL2',data_file_level2,'QNP',4096,istat) !opens rz if(istat.ne.0) goto 19 call book_level2 if(DEBUG)then print*,'(creating DEBUG nt-uple and histos)' call book_debug endif c------------------------------------------------------------------------ c c fills bad variables with online DSP bad strips from calib histograms c for bad strip exclusion c c------------------------------------------------------------------------ # ifdef TEST2003 c------------------------------------------------------------------------ c c opens calib file c c------------------------------------------------------------------------ 505 format(a,'output_',a,'_calib.rz') write(data_file_calib,505) $ data_dir(1:LNBLNK(data_dir)) $ ,data_file(1:LNBLNK(data_file)) print*,' ' print*,'OPENING CALIB FILE...', data_file_calib IQUEST(10)=65000 call HROPEN(lun_data_calib,'IN',data_file_calib,'QP',4096,istat) if(istat.ne.0) goto 17 call HCDIR('//IN',' ') call HRIN(0,9999,0) !puts histograms in memory print*,' ' print*,' ' print*,'READING PEDESTAL, SIGMA AND BADSTRIP HISTOGRAMS...' print*,' ' call fillpedsig call HREND('IN') close (lun_data_calib) # else 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 # endif c------------------------------------------------------------------------ c c opens level1 file c c------------------------------------------------------------------------ # ifdef TEST2003 504 format(a,'output_',a,'_level1.rz') # else 504 format(a,'DW_',a,'_level1.rz') # endif write(data_file_level1,504) $ data_dir(1:LNBLNK(data_dir)) $ ,data_file(1:LNBLNK(data_file)) print*,'' print*,'OPENING LEVEL1 FILE:' print*,data_file_level1 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,'QP',4096,istat) !opens rz if(istat.ne.0) goto 19 call HRIN(ntp_level1,9999,20) call access_level1 c call HPRNTU(ntp_level1+20) call HNOENT(ntp_level1+20,iemax0) print*,' events',iemax0 ************************************************************ ************************************************************ ************************************************************ * * start track analysis * ************************************************************ ************************************************************ ************************************************************ maxevent=minevent+ntotev-1 do iev = minevent,MIN(iemax0,maxevent) !loop on events call HCDIR('//LEVEL1',' ') call HGNT(ntp_level1+20,iev,ierr) !reads an event if(ierr.ne.0) goto 21 *------------------------------------------------------ * LEVEL2 N-TUPLE INITIALIZATIONS call init_level2 if(DEBUG)call init_level2_debug if(DEBUG)DBG_FILLED=.false. if(.not.good1)then goto 8800 !fill nt-uple and go to next event endif *------------------------------------------------------ # ifndef TEST2003 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 ',nev2 $ ,' (CPU pkt N.',pkt_num1,')' print*,'--> ',data_file_calib IQUEST(10)=65000 call HROPEN(lun_data_calib, $ '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_calib) which_calib_last=which_calib elseif(which_calib.eq.0)then nocalib=nocalib+1 good2=.false. goto 8800 !fill nt-uple and go to next event endif # endif *------------------------------------------------------ * cut on maximum number of clusters *------------------------------------------------------ if(nclstr1.gt.nclstrmax_level2)then goto 8800 !fill nt-uple and go to next event endif do i=1,nclstr1 cl_used(i)=0 !init mask of clusters associated to a track enddo if(DEBUG)then print*,'----------------------------------' print*,iev,' ** ',nev2 endif *------------------------------------------------------------------------------- * STEP 1 *------------------------------------------------------------------------------- * X-Y cluster association * * Clusters are associated to form COUPLES * Clusters not associated in any couple are called SINGLETS * * Track identification (Hough transform) and fitting is first done on couples. * Hence singlets are possibly added to the track. * * Variables assigned by the routine "cl_to_couples" are those in the * common blocks: * - common/clusters/cl_good * - common/couples/clx,cly,ncp_plane,ncp_tot,cp_useds1,cp_useds2 * - common/singlets/ncls,cls,cl_single *------------------------------------------------------------------------------- *------------------------------------------------------------------------------- iflag=0 call cl_to_couples(iflag) if(iflag.eq.1)then !bad event goto 880 !fill ntp and go to next event endif *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * selezione di tracce pulite per diagnostica *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c$$$ if(DEBUG)then c$$$ do ip=1,nplanes c$$$ if(ncp_plane(ip).ne.1)good2=.false. c$$$ enddo c$$$c if(good2.eq.0)goto 100!next event c$$$c if(good2.eq.0)goto 880!fill ntp and go to next event c$$$ endif *----------------------------------------------------- *----------------------------------------------------- * HOUGH TRASFORM *----------------------------------------------------- *----------------------------------------------------- *------------------------------------------------------------------------------- * STEP 2 *------------------------------------------------------------------------------- * * Association of couples to form * - DOUBLETS in YZ view * - TRIPLETS in XZ view * * Variables assigned by the routine "cp_to_doubtrip" are those in the * common blocks: * - common/hough_param/ * $ alfayz1, !Y0 * $ alfayz2, !tg theta-yz * $ alfaxz1, !X0 * $ alfaxz2, !tg theta-xz * $ alfaxz3 !1/r * - common/doublets/ndblt,cpyz1,cpyz2 * - common/triplets/ntrpt,cpxz1,cpxz2,cpxz3 *------------------------------------------------------------------------------- *------------------------------------------------------------------------------- iflag=0 call cp_to_doubtrip(iflag) if(iflag.eq.1)then !bad event goto 880 !fill ntp and go to next event endif *------------------------------------------------------------------------------- * STEP 3 *------------------------------------------------------------------------------- * * Classification of doublets and triplets to form CLOUDS, * according to distance in parameter space. * * cloud = cluster of points (doublets/triplets) in parameter space * * * * Variables assigned by the routine "doub_to_YZcloud" are those in the * common blocks: * - common/clouds_yz/ * $ nclouds_yz * $ ,alfayz1_av,alfayz2_av * $ ,ptcloud_yz,db_cloud,cpcloud_yz * * Variables assigned by the routine "trip_to_XZcloud" are those in the * common blocks: * common/clouds_xz/ * $ nclouds_xz xz2_av,alfaxz3_av * $ ,ptcloud_xz,tr_cloud,cpcloud_xz *------------------------------------------------------------------------------- *------------------------------------------------------------------------------- iflag=0 call doub_to_YZcloud(iflag) if(iflag.eq.1)then !bad event goto 880 !fill ntp and go to next event endif iflag=0 call trip_to_XZcloud(iflag) if(iflag.eq.1)then !bad event goto 880 !fill ntp and go to next event endif c***************************************************** cccccc 01/12/2005 modified by elena if(DEBUG)then call fill_level2_clouds call HCDIR('//LEVEL2',' ') call HFNT(ntp_level2+1) !fill DEBUG nt-uple DBG_FILLED=.true. endif c***************************************************** *------------------------------------------------------------------------------- * STEP 4 (ITERATED until any other physical track isn't found) *------------------------------------------------------------------------------- * * YZ and XZ clouds are combined in order to obtain the initial guess * of the candidate-track parameters. * A minimum number of matching couples between YZ and XZ clouds is required. * * A TRACK CANDIDATE is defined by * - the couples resulting from the INTERSECTION of the two clouds, and * - the associated track parameters (evaluated by performing a zero-order * track fitting) * * The NTRACKS candidate-track parameters are stored in common block: * * - common/track_candidates/NTRACKS,AL_STORE * $ ,XV_STORE,YV_STORE,ZV_STORE * $ ,XM_STORE,YM_STORE,ZM_STORE * $ ,RESX_STORE,RESY_STORE * $ ,AXV_STORE,AYV_STORE * $ ,XGOOD_STORE,YGOOD_STORE * $ ,CP_STORE,RCHI2_STORE * *------------------------------------------------------------------------------- *------------------------------------------------------------------------------- ntrk=0 !counter of identified physical tracks 11111 continue !<<<<<<< come here when performing a new search iflag=0 call clouds_to_ctrack(iflag) if(iflag.eq.1)then !no candidate tracks found goto 880 !fill ntp and go to next event endif FIMAGE=.false. !processing best track (not track image) ibest=0 !best track among candidates iimage=0 !track image * ------------- select the best track ------------- rchi2best=1000000000. do i=1,ntracks if(RCHI2_STORE(i).lt.rchi2best.and. $ RCHI2_STORE(i).gt.0)then ibest=i rchi2best=RCHI2_STORE(i) endif enddo if(ibest.eq.0)goto 880 !>> no good candidates *------------------------------------------------------------------------------- * The best track candidate (ibest) is selected and a new fitting is performed. * Previous to this, the track is refined by: * - possibly adding new COUPLES or SINGLETS from the missing planes * - evaluating the coordinates with improved PFAs * ( angle-dependent ETA algorithms ) *------------------------------------------------------------------------------- 1212 continue !<<<<< come here to fit track-image if(.not.FIMAGE)then !processing best candidate icand=ibest else !processing image icand=iimage iimage=0 endif if(icand.eq.0)then print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand $ ,ibest,iimage return endif * *-*-*-*-*-*-*-*-*-*-*-*-*-*-* call refine_track(icand) * *-*-*-*-*-*-*-*-*-*-*-*-*-*-* * ********************************************************** * ************************** FIT *** FIT *** FIT *** FIT *** * ********************************************************** do i=1,5 AL(i)=dble(AL_STORE(i,icand)) enddo ifail=0 !error flag in chi2 computation jstep=0 !# minimization steps call mini_2(jstep,ifail) if(ifail.ne.0) then if(DEBUG)then print *, $ '*** MINIMIZATION FAILURE *** (mini_2) ' $ ,iev endif chi2=-chi2 endif if(DEBUG)then print*,'----------------------------- improved track coord' 22222 format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5) do ip=1,6 write(*,22222)ip,zm(ip),xm(ip),ym(ip) $ ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip) $ ,xgood(ip),ygood(ip),resx(ip),resy(ip) enddo endif c rchi2=chi2/dble(ndof) if(DEBUG)then print*,' ' print*,'****** SELECTED TRACK *************' print*,'# R. chi2 RIG' print*,' --- ',chi2,' --- ' $ ,1./abs(AL(5)) print*,'***********************************' endif * ********************************************************** * ************************** FIT *** FIT *** FIT *** FIT *** * ********************************************************** * ------------- search if the track has an IMAGE ------------- * ------------- (also this is stored ) ------------- if(FIMAGE)goto 122 !>>> jump! (this is already an image) * now search for track-image, by comparing couples IDs do i=1,ntracks iimage=i do ip=1,nplanes if( CP_STORE(nplanes-ip+1,icand).ne. $ -1*CP_STORE(nplanes-ip+1,i) )iimage=0 enddo if( iimage.ne.0.and. c $ RCHI2_STORE(i).le.CHI2MAX.and. c $ RCHI2_STORE(i).gt.0.and. $ .true.)then if(DEBUG)print*,'Track candidate ',iimage $ ,' >>> TRACK IMAGE >>> of' $ ,ibest goto 122 !image track found endif enddo 122 continue * --- and store the results -------------------------------- ntrk = ntrk + 1 !counter of found tracks if(.not.FIMAGE $ .and.iimage.eq.0) image(ntrk)= 0 if(.not.FIMAGE $ .and.iimage.ne.0)image(ntrk)=ntrk+1 !this is the image of the next if(FIMAGE) image(ntrk)=ntrk-1 !this is the image of the previous call fill_level2_tracks(ntrk) !==> good2=.true. c print*,'++++++++++ iimage,fimage,ntrk,image ' c $ ,iimage,fimage,ntrk,image(ntrk) if(ntrk.eq.NTRKMAX)then if(DEBUG) $ print*, $ '** warning ** number of identified '// $ 'tracks exceeds vector dimension ' $ ,'( ',NTRKMAX,' )' cc good2=.false. goto 880 !fill ntp and go to next event endif if(iimage.ne.0)then FIMAGE=.true. ! goto 1212 !>>> fit image-track endif * --- then remove selected clusters (ibest+iimage) from clouds ---- call clean_XYclouds(ibest) if(iflag.eq.1)then !bad event goto 880 !fill ntp and go to next event endif * ********************************************************** * condition to start a new search * ********************************************************** ixznew=0 do ixz=1,nclouds_xz if(ptcloud_xz(ixz).ge.nptxz_min)ixznew=1 enddo iyznew=0 do iyz=1,nclouds_yz if(ptcloud_yz(iyz).ge.nptyz_min)iyznew=1 enddo if(ixznew.ne.0.and. $ iyznew.ne.0.and. $ rchi2best.le.CHI2MAX.and. c $ rchi2best.lt.15..and. $ .true.)then if(DEBUG)then print*,'***** NEW SEARCH ****' endif goto 11111 !try new search endif * ********************************************** * >>>>>>>>>>>>>>>>>>> NT-UPLE filling <<<<<<<<<<<<<<<<<<<< * + + + + + + + + + + + + + + + + + * / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / * + + + + + + + + + + + + + + + + + 880 continue * ********************************************************** * stores info about clusters not associated with any track * ********************************************************** c***************************************************** cccccc 27/9/2005 modified by david fedele * count #cluster per plane not associated to any track c$$$ do icl=1,nclstr1 c$$$ if(cl_used(icl).eq.0)then !cluster not included in any track c$$$ ip=nplanes-npl(VIEW(icl))+1 c$$$ if(mod(VIEW(icl),2).eq.0)nclsx(ip)=nclsx(ip)+1 c$$$ if(mod(VIEW(icl),2).eq.1)nclsy(ip)=nclsy(ip)+1 c$$$ endif c$$$c print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip) c$$$ enddo c********************************************************************** call fill_level2_siglets c print*,'****** ',iev,' - ',ntrk,nclsx,nclsy c print*,'****** ',iev,' - ',ntrk,' --- ',(BdL(i),i=1,ntrk) if(DEBUG)then print*,'' print*,'DONE!' print*,'' print*,'* summary *' print*,'tracks ',ntrk print*,'cl used ',(cl_used(i),i=1,nclstr1) c***************************************************** c$$$cccccc 27/9/2005 modified by david fedele c$$$ print*,'cl unused (x-y)' c$$$ do ip=1,nplanes c$$$ print*,ip,' << ',nclsx(ip),nclsy(ip) c$$$ enddo c****************************************************** print*,'' print*,'' endif 8800 continue call HCDIR('//LEVEL2',' ') call HFNT(ntp_level2) !fill LEVEL2 nt-uple if(.not.DBG_FILLED.and.DEBUG) $ print*,'@@@@@@@@@@@@',ntp_level2+1,nev2_nt if(.not.DBG_FILLED.and.DEBUG)call HFNT(ntp_level2+1) 100 continue enddo !end loop on events c------------------------------------------------------------------------ c c no error exit c c------------------------------------------------------------------------ c$$$ print*,' ' c$$$ print*,'REDUCTION SUCCESSFULLY COMPLETED' c$$$ print*,' ' c$$$ print*,' ' goto 9000 !happy ending c------------------------------------------------------------------------ c c data file opening error c c------------------------------------------------------------------------ 19 continue print*,' ' print*,'ERROR OPENING DATA FILE: ',data_file print*,' ' print*,' ' goto 9000 !the end c------------------------------------------------------------------------ c c level1 ntuple event reading error c c------------------------------------------------------------------------ 21 continue print*,' ' print*,'ERROR WHILE READING LEVEL1 NTUPLE, AT EVENT $ : ',iev print*,' ' print*,' ' goto 9000 !the end c------------------------------------------------------------------------ c c calib file opening error c c------------------------------------------------------------------------ 17 continue print*,' ' print*,'preanalysis: ERROR OPENING INPUT FILE: ',data_file_calib print*,' ' print*,' ' goto 9000 !the end c------------------------------------------------------------------------ c c closes files and exits c c------------------------------------------------------------------------ 9000 continue if(DEBUG)call HPRNTU(ntp_level2+1) call HPRNTU(ntp_level2) print*,'' call HCDIR('//LEVEL2',' ') call HROUT(ntp_level2,ICYCLE,'T') print *,'- Stored LEVEL2 nt-uple (',nev2,' entries )' if(DEBUG)call HROUT(ntp_level2+1,ICYCLE,'T') if(DEBUG)print *,'- Stored DEBUG nt-uple ' call HREND('level2') close(lun_data_level2) call HCDIR('//LEVEL1',' ') call HREND('level1') close(lun_data_level1) stop end ************************************************************ # include "momanhough-subroutines.F"