************************************************************************* * 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' include '../common/momanhough_init.f' 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' * flag to set debug mode logical DEBUG common/dbg/DEBUG * 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 ! * ****************************************** * PLANE TO EXCLUDE (for efficiency studies) * --> provided as input parameter * (set according to "tracking" notation): * 1-6 from TOP to BOTTOM * ****************************************** integer ipexclude common/plane_exclude/ipexclude COMMON/QUEST/IQUEST(100) c !permette di ottenere ntuple funzionanti nonostante c ! il messaggio dei 64K di RZOUT...!??? 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*,'Plane to exclude from fitting (1-6 from TOP to BOTTOM):' read(*,*) ipexclude print*, ipexclude print*,'---------------------------------------------------' ****** INITIALIZATIONS ************************************* * 1) read charge-correlation parameters print*,' ' print*,'- read charge-correlation parameters' print*,' ' call readchargeparam * 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 PFA parameters' print*,' ' call readetaparam * 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_no-p',i1,'.rz') write(data_file_level2,503) $ data_dir(1:LNBLNK(data_dir)) $ ,data_file(1:LNBLNK(data_file)) $ ,ipexclude 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(.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 * $ ,alfaxz1_av,alfaxz2_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 *------------------------------------------------------------------------------- * 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=100000. 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 track 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) * *************************************** * Checks if a plane should be excluded * (for efficiency study) * *************************************** xgood(ipexclude) = 0 !<<< ygood(ipexclude) = 0 !<<< * *************************************** * ********************************************************** * ************************** 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,' * ',3f9.4,' --- ',4f9.4,' --- ',2f4.0,2f8.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 * count #cluster per plane not associated to any track do icl=1,nclstr1 if(cl_used(icl).eq.0)then !cluster not included in any track ip=nplanes-npl(VIEW(icl))+1 if(mod(VIEW(icl),2).eq.0)nclsx(ip)=nclsx(ip)+1 if(mod(VIEW(icl),2).eq.1)nclsy(ip)=nclsy(ip)+1 endif c print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip) enddo if(DEBUG)then print*,'' print*,'DONE!' print*,'' print*,'* summary *' print*,'tracks ',ntrk print*,'cl used ',(cl_used(i),i=1,nclstr1) print*,'cl unused (x-y)' do ip=1,nplanes print*,ip,' << ',nclsx(ip),nclsy(ip) enddo print*,'' print*,'' endif 8800 continue call HCDIR('//LEVEL2',' ') call HFNT(ntp_level2) !fill LEVEL2 nt-uple if(DEBUG)then good2_nt=good2 nev2_nt=nev2 if(ntrpt.gt.ntrpt_max_nt.or.ndblt.gt.ndblt_max_nt)then good2_nt=.false. ntrpt_nt=0 ndblt_nt=0 NCLOUDS_XZ=0 NCLOUDS_YZ=0 else ndblt_nt=ndblt ntrpt_nt=ntrpt do id=1,ndblt alfayz1_nt(id)=alfayz1(id) !Y0 alfayz2_nt(id)=alfayz2(id) !tg theta-yz db_cloud_nt(id)=db_cloud(id) enddo do it=1,ntrpt alfaxz1_nt(it)=alfaxz1(it) !X0 alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz alfaxz3_nt(it)=alfaxz3(it) !1/r tr_cloud_nt(it)=tr_cloud(it) enddo nclouds_yz_nt=nclouds_yz nclouds_xz_nt=nclouds_xz do iyz=1,nclouds_yz ptcloud_yz_nt(iyz)=ptcloud_yz(iyz) alfayz1_av_nt(iyz)=alfayz1_av(iyz) alfayz2_av_nt(iyz)=alfayz2_av(iyz) enddo do ixz=1,nclouds_xz ptcloud_xz_nt(ixz)=ptcloud_xz(ixz) alfaxz1_av_nt(ixz)=alfaxz1_av(ixz) alfaxz2_av_nt(ixz)=alfaxz2_av(ixz) alfaxz3_av_nt(ixz)=alfaxz3_av(ixz) enddo endif c print*,ntrpt,ntrpt_max_nt,ndblt,ndblt_max_nt c do iii=1,ndblt c print*,iii,alfayz1(iii),alfayz2(iii) c enddo c print*,good2_nt,good2 call HFNT(ntp_level2+1) !fill DEBUG nt-uple endif 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" c$$$ c$$$ c$$$ c$$$************************************************************ c$$$ c$$$ subroutine readchargeparam c$$$ c$$$ c$$$ include '../common/commontracker.f' c$$$ include '../common/calib.f' c$$$ c$$$c character*40 fname_binning c$$$ character*60 fname_param c$$$c character*120 cmd1 c$$$c character*120 cmd2 c$$$ c$$$c 201 format('../common/charge_param/charge-l',i1,'.dat') c$$$c 201 format('$TRK_GRND/source/common/charge_param/charge-l',i1,'.dat') c$$$ 201 format('charge-l',i1,'.dat') c$$$ do ilad=1,nladders_view c$$$ write(fname_param,201)ilad c$$$c$$$ cmd1='cp $TRK_GRND/source/common/charge_param/' c$$$c$$$ $ //fname_param(1:LNBLNK(fname_param))//' .' c$$$c$$$ call system(cmd1) c$$$ print *,'Opening file: ',fname_param c$$$ open(10, c$$$ $ FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param)) c$$$ $ ,STATUS='UNKNOWN' c$$$ $ ,IOSTAT=iostat c$$$ $ ) c$$$ if(iostat.ne.0)then c$$$ print*,'READCHARGEPARAM: *** Error in opening file ***' c$$$ return c$$$ endif c$$$ do ip=1,nplanes c$$$ read(10,* c$$$ $ ,IOSTAT=iostat c$$$ $ )pip, c$$$ $ kch(ip,ilad),cch(ip,ilad),sch(ip,ilad) c$$$c print*,ilad,ip,pip,kch(ip,ilad), c$$$c $ cch(ip,ilad),sch(ip,ilad) c$$$ enddo c$$$ close(10) c$$$c$$$ cmd2='rm -f ' c$$$c$$$ $ //fname_param(1:LNBLNK(fname_param)) c$$$c$$$ call system(cmd2) c$$$ enddo c$$$ c$$$ return c$$$ end c$$$ c$$$************************************************************ c$$$************************************************************ c$$$************************************************************ c$$$************************************************************ c$$$* c$$$* This routine provides the coordinates (in cm) in the PAMELA reference system: c$$$* - of the point associated with a COUPLE ---> (xPAM,yPAM,zPAM) c$$$* - of the extremes of the segment c$$$* associated with a SINGLET ---------------> (xPAM_A,yPAM_A,zPAM_A) c$$$* ---> (xPAM_B,yPAM_B,zPAM_B) c$$$* c$$$* It also assigns the spatial resolution to the evaluated coordinates, c$$$* as a function (in principle) of the multiplicity, the angle, the PFA etc... c$$$* c$$$* c$$$* To call the routine you must pass the arguments: c$$$* icx - ID of cluster x c$$$* icy - ID of cluster y c$$$* sensor - sensor (1,2) c$$$* PFAx - Position Finding Algorithm in x (COG2,ETA2,...) c$$$* PFAy - Position Finding Algorithm in y (COG2,ETA2,...) c$$$* angx - Projected angle in x c$$$* angy - Projected angle in y c$$$* c$$$* --------- COUPLES ------------------------------------------------------- c$$$* The couple defines a point in the space. c$$$* The coordinates of the point are evaluated as follows: c$$$* 1 - the corrected coordinates relative to the sensor are evaluated c$$$* according to the chosen PFA --> (xi,yi,0) c$$$* 2 - coordinates are rotated and traslated, according to the aligmnet c$$$* parameters, and expressed in the reference system of the mechanical c$$$* sensor --> (xrt,yrt,zrt) c$$$* 3 - coordinates are finally converted to the PAMELA reference system c$$$* --> (xPAM,yPAM,zPAM) c$$$* c$$$* --------- SINGLETS ------------------------------------------------------- c$$$* Since a coordinate is missing, the singlet defines not a point c$$$* in the space but a segment AB (parallel to the strips). c$$$* In this case the routine returns the coordinates in the PAMELA reference c$$$* system of the two extremes A and B of the segment: c$$$* --> (xPAM_A,yPAM_A,zPAM_A) c$$$* --> (xPAM_B,yPAM_B,zPAM_B) c$$$* c$$$* ========================================================== c$$$* c$$$* The output of the routine is stored in the commons: c$$$* c$$$* double precision xPAM,yPAM,zPAM c$$$* common/coord_xyz_PAM/xPAM,yPAM,zPAM c$$$* c$$$* double precision xPAM_A,yPAM_A,zPAM_A c$$$* double precision xPAM_B,yPAM_B,zPAM_B c$$$* common/coord_AB_PAM/xPAM_A,yPAM_A,zPAM_A,xPAM_B,yPAM_B,zPAM_B c$$$* c$$$* double precision resxPAM,resyPAM c$$$* common/resolution_PAM/resxPAM,resyPAM c$$$* c$$$* (in file ../common/common_xyzPAM.f) c$$$* c$$$* c$$$ c$$$ subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,angx,angy) c$$$ c$$$ include '../common/commontracker.f' c$$$ include '../common/calib.f' c$$$ include '../common/level1.f' c$$$ include '../common/common_align.f' c$$$ include '../common/common_mech.f' c$$$ include '../common/common_xyzPAM.f' c$$$ include '../common/common_resxy.f' c$$$ c$$$ logical DEBUG c$$$ common/dbg/DEBUG c$$$ c$$$ integer icx,icy !X-Y cluster ID c$$$ integer sensor c$$$ integer viewx,viewy c$$$ character*4 PFAx,PFAy !PFA to be used c$$$ real angx,angy !X-Y angle c$$$ c$$$ real stripx,stripy c$$$ c$$$ double precision xrt,yrt,zrt c$$$ double precision xrt_A,yrt_A,zrt_A c$$$ double precision xrt_B,yrt_B,zrt_B c$$$c double precision xi,yi,zi c$$$c double precision xi_A,yi_A,zi_A c$$$c double precision xi_B,yi_B,zi_B c$$$ c$$$ c$$$ parameter (ndivx=30) c$$$ c$$$ c$$$ resxPAM = 0 c$$$ resyPAM = 0 c$$$ c$$$ xPAM = 0. c$$$ yPAM = 0. c$$$ zPAM = 0. c$$$ xPAM_A = 0. c$$$ yPAM_A = 0. c$$$ zPAM_A = 0. c$$$ xPAM_B = 0. c$$$ yPAM_B = 0. c$$$ zPAM_B = 0. c$$$ c$$$* ----------------- c$$$* CLUSTER X c$$$* ----------------- c$$$ c$$$ if(icx.ne.0)then c$$$ viewx = VIEW(icx) c$$$ nldx = nld(MAXS(icx),VIEW(icx)) c$$$ nplx = npl(VIEW(icx)) c$$$ resxPAM = RESXAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!! c$$$ c$$$ stripx = float(MAXS(icx)) c$$$ if(PFAx.eq.'COG2')then c$$$ stripx = stripx + cog(2,icx) c$$$ resxPAM = resxPAM*fbad_cog(2,icx) c$$$ elseif(PFAx.eq.'ETA2')then c$$$ cog2 = cog(2,icx) c$$$ etacorr = pfa_eta2(cog2,viewx,nldx,angx) c$$$c print*,'---- ',viewx,cog2,etacorr c$$$ stripx = stripx + etacorr c$$$c print*,'COG2 ',cog2,' >>> ETA2 ',etacorr c$$$ if(DEBUG.and.fbad_cog(2,icx).ne.1) c$$$ $ print*,'BAD icx >>> ',viewx,fbad_cog(2,icx) c$$$ resxPAM = resxPAM*fbad_cog(2,icx) c$$$ else c$$$ print*,'Accontentati di MAXS !!!!' c$$$ endif c$$$ c$$$ c$$$ endif c$$$ c$$$* ----------------- c$$$* CLUSTER Y c$$$* ----------------- c$$$ c$$$ if(icy.ne.0)then c$$$ viewy = VIEW(icy) c$$$ nldy = nld(MAXS(icy),VIEW(icy)) c$$$ nply = npl(VIEW(icy)) c$$$ resyPAM = RESYAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!! c$$$ c$$$ c$$$ if(icx.ne.0.and.(nply.ne.nplx.or.nldy.ne.nldx))then c$$$ print*,'xyz_PAM ***ERROR*** invalid cluster couple!!! ' c$$$ $ ,icx,icy c$$$ goto 100 c$$$ endif c$$$ c$$$ stripy = float(MAXS(icy)) c$$$ if(PFAy.eq.'COG2')then c$$$ stripy = stripy + cog(2,icy) c$$$ resyPAM = resyPAM*fbad_cog(2,icy) c$$$ elseif(PFAy.eq.'ETA2')then c$$$ cog2 = cog(2,icy) c$$$ etacorr = pfa_eta2(cog2,viewy,nldy,angy) c$$$ stripy = stripy + etacorr c$$$c print*,'COG2 ',cog2,' >>> ETA2 ',etacorr c$$$ resyPAM = resyPAM*fbad_cog(2,icy) c$$$ if(DEBUG.and.fbad_cog(2,icy).ne.1) c$$$ $ print*,'BAD icy >>> ',viewy,fbad_cog(2,icy) c$$$ else c$$$ print*,'Accontentati di MAXS !!!!' c$$$ endif c$$$ c$$$ endif c$$$ c$$$ c$$$c=========================================================== c$$$C COUPLE c$$$C=========================================================== c$$$ if(icx.ne.0.and.icy.ne.0)then c$$$ c$$$c------------------------------------------------------------------------ c$$$c (xi,yi,zi) = mechanical coordinates in the silicon sensor frame c$$$c------------------------------------------------------------------------ c$$$ xi = acoordsi(stripx,viewx) c$$$ yi = acoordsi(stripy,viewy) c$$$ zi = 0. c$$$ c$$$ c$$$c------------------------------------------------------------------------ c$$$c (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame c$$$c------------------------------------------------------------------------ c$$$c N.B. I convert angles from microradiants to radiants c$$$ c$$$ xrt = xi c$$$ $ - omega(nplx,nldx,sensor)*yi c$$$ $ + gamma(nplx,nldx,sensor)*zi c$$$ $ + dx(nplx,nldx,sensor) c$$$ c$$$ yrt = omega(nplx,nldx,sensor)*xi c$$$ $ + yi c$$$ $ - beta(nplx,nldx,sensor)*zi c$$$ $ + dy(nplx,nldx,sensor) c$$$ c$$$ zrt = -gamma(nplx,nldx,sensor)*xi c$$$ $ + beta(nplx,nldx,sensor)*yi c$$$ $ + zi c$$$ $ + dz(nplx,nldx,sensor) c$$$ c$$$c xrt = xi c$$$c yrt = yi c$$$c zrt = zi c$$$ c$$$c------------------------------------------------------------------------ c$$$c (xPAM,yPAM,zPAM) = measured coordinates (in cm) c$$$c in PAMELA reference system c$$$c------------------------------------------------------------------------ c$$$ c$$$ xPAM = dcoord(xrt,viewx,nldx,sensor) / 1.d4 c$$$ yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4 c$$$ zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4 c$$$ c$$$ xPAM_A = 0. c$$$ yPAM_A = 0. c$$$ zPAM_A = 0. c$$$ c$$$ xPAM_B = 0. c$$$ yPAM_B = 0. c$$$ zPAM_B = 0. c$$$ c$$$ elseif( c$$$ $ (icx.ne.0.and.icy.eq.0).or. c$$$ $ (icx.eq.0.and.icy.ne.0).or. c$$$ $ .false. c$$$ $ )then c$$$ c$$$c------------------------------------------------------------------------ c$$$c (xi,yi,zi) = mechanical coordinates in the silicon sensor frame c$$$c------------------------------------------------------------------------ c$$$ c$$$ if(icy.ne.0)then c$$$c=========================================================== c$$$C Y-SINGLET c$$$C=========================================================== c$$$ nplx = nply c$$$ nldx = nldy c$$$ viewx = viewy + 1 c$$$ c$$$ yi = acoordsi(stripy,viewy) c$$$ c$$$ xi_A = edgeY_d - SiDimX/2 c$$$ yi_A = yi c$$$ zi_A = 0. c$$$ c$$$ xi_B = SiDimX/2 - edgeY_u c$$$ yi_B = yi c$$$ zi_B = 0. c$$$ c$$$c print*,'Y-cl ',icy,stripy,' --> ',yi c$$$c print*,xi_A,' <--> ',xi_B c$$$ c$$$ elseif(icx.ne.0)then c$$$c=========================================================== c$$$C X-SINGLET c$$$C=========================================================== c$$$ c$$$ nply = nplx c$$$ nldy = nldx c$$$ viewy = viewx - 1 c$$$ c$$$ xi = acoordsi(stripx,viewx) c$$$ c$$$ xi_A = xi c$$$ yi_A = edgeX_d - SiDimY/2 c$$$ zi_A = 0. c$$$ c$$$ xi_B = xi c$$$ yi_B = SiDimY/2 - edgeX_u c$$$ zi_B = 0. c$$$ c$$$ if(viewy.eq.11)then c$$$ yi = yi_A c$$$ yi_A = yi_B c$$$ yi_B = yi c$$$ endif c$$$ c$$$c print*,'X-cl ',icx,stripx,' --> ',xi c$$$c print*,yi_A,' <--> ',yi_B c$$$ c$$$ else c$$$ c$$$ print *,'routine xyz_PAM ---> not properly used !!!' c$$$ print *,'icx = ',icx c$$$ print *,'icy = ',icy c$$$ goto 100 c$$$ c$$$ endif c$$$c------------------------------------------------------------------------ c$$$c (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame c$$$c------------------------------------------------------------------------ c$$$c N.B. I convert angles from microradiants to radiants c$$$ c$$$ xrt_A = xi_A c$$$ $ - omega(nplx,nldx,sensor)*yi_A c$$$ $ + gamma(nplx,nldx,sensor)*zi_A c$$$ $ + dx(nplx,nldx,sensor) c$$$ c$$$ yrt_A = omega(nplx,nldx,sensor)*xi_A c$$$ $ + yi_A c$$$ $ - beta(nplx,nldx,sensor)*zi_A c$$$ $ + dy(nplx,nldx,sensor) c$$$ c$$$ zrt_A = -gamma(nplx,nldx,sensor)*xi_A c$$$ $ + beta(nplx,nldx,sensor)*yi_A c$$$ $ + zi_A c$$$ $ + dz(nplx,nldx,sensor) c$$$ c$$$ xrt_B = xi_B c$$$ $ - omega(nplx,nldx,sensor)*yi_B c$$$ $ + gamma(nplx,nldx,sensor)*zi_B c$$$ $ + dx(nplx,nldx,sensor) c$$$ c$$$ yrt_B = omega(nplx,nldx,sensor)*xi_B c$$$ $ + yi_B c$$$ $ - beta(nplx,nldx,sensor)*zi_B c$$$ $ + dy(nplx,nldx,sensor) c$$$ c$$$ zrt_B = -gamma(nplx,nldx,sensor)*xi_B c$$$ $ + beta(nplx,nldx,sensor)*yi_B c$$$ $ + zi_B c$$$ $ + dz(nplx,nldx,sensor) c$$$ c$$$ c$$$c xrt = xi c$$$c yrt = yi c$$$c zrt = zi c$$$ c$$$c------------------------------------------------------------------------ c$$$c (xPAM,yPAM,zPAM) = measured coordinates (in cm) c$$$c in PAMELA reference system c$$$c------------------------------------------------------------------------ c$$$ c$$$ xPAM = 0. c$$$ yPAM = 0. c$$$ zPAM = 0. c$$$ c$$$ xPAM_A = dcoord(xrt_A,viewx,nldx,sensor) / 1.d4 c$$$ yPAM_A = dcoord(yrt_A,viewy,nldy,sensor) / 1.d4 c$$$ zPAM_A = ( zrt_A + z_mech_sensor(nplx,nldx,sensor)*1000.)/ 1.d4 c$$$ c$$$ xPAM_B = dcoord(xrt_B,viewx,nldx,sensor) / 1.d4 c$$$ yPAM_B = dcoord(yrt_B,viewy,nldy,sensor) / 1.d4 c$$$ zPAM_B = ( zrt_B + z_mech_sensor(nplx,nldx,sensor)*1000.)/ 1.d4 c$$$ c$$$ c$$$c print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')' c$$$ c$$$ else c$$$ c$$$ print *,'routine xyz_PAM ---> not properly used !!!' c$$$ print *,'icx = ',icx c$$$ print *,'icy = ',icy c$$$ c$$$ endif c$$$ c$$$ 100 continue c$$$ end c$$$ c$$$ c$$$******************************************************************************** c$$$******************************************************************************** c$$$******************************************************************************** c$$$* c$$$* The function distance_to(XP,YP) should be used after c$$$* a call to the xyz_PAM routine and it evaluate the c$$$* NORMALIZED distance (PROJECTED on the XY plane) between c$$$* the point (XP,YP), argument of the function, c$$$* and: c$$$* c$$$* - the point (xPAM,yPAM,zPAM), in the case of a COUPLE c$$$* or c$$$* - the segment (xPAM_A,yPAM_A,zPAM_A)-(xPAM_B,yPAM_B,zPAM_B), c$$$* in the case of a SINGLET. c$$$* c$$$* ( The routine xyz_PAM fills the common defined in "common_xyzPAM.f", c$$$* which stores the coordinates of the couple/singlet ) c$$$* c$$$******************************************************************************** c$$$ c$$$ real function distance_to(XPP,YPP) c$$$ c$$$ include '../common/common_xyzPAM.f' c$$$ c$$$* ----------------------------------- c$$$* it computes the normalized distance c$$$* ( i.e. distance/resolution ) c$$$* ----------------------------------- c$$$ c$$$ double precision distance,RE c$$$ double precision BETA,ALFA,xmi,ymi c$$$ c$$$* ---------------------- c$$$ if ( c$$$ + xPAM.eq.0.and. c$$$ + yPAM.eq.0.and. c$$$ + zPAM.eq.0.and. c$$$ + xPAM_A.ne.0.and. c$$$ + yPAM_A.ne.0.and. c$$$ + zPAM_A.ne.0.and. c$$$ + xPAM_B.ne.0.and. c$$$ + yPAM_B.ne.0.and. c$$$ + zPAM_B.ne.0.and. c$$$ + .true.)then c$$$* ----------------------- c$$$* DISTANCE TO --- SINGLET c$$$* ----------------------- c$$$ if(abs(sngl(xPAM_B-xPAM_A)).lt.abs(sngl(yPAM_B-yPAM_A)))then c$$$* |||---------- X CLUSTER c$$$ c$$$ BETA = (xPAM_B-xPAM_A)/(yPAM_B-yPAM_A) c$$$ ALFA = xPAM_A - BETA * yPAM_A c$$$ c$$$ ymi = ( YPP + BETA*XPP - BETA*ALFA )/(1+BETA**2) c$$$ if(ymi.lt.dmin1(yPAM_A,yPAM_B))ymi=dmin1(yPAM_A,yPAM_B) c$$$ if(ymi.gt.dmax1(yPAM_A,yPAM_B))ymi=dmax1(yPAM_A,yPAM_B) c$$$ xmi = ALFA + BETA * ymi c$$$ RE = resxPAM c$$$ c$$$ else c$$$* |||---------- Y CLUSTER c$$$ c$$$ BETA = (yPAM_B-yPAM_A)/(xPAM_B-xPAM_A) c$$$ ALFA = yPAM_A - BETA * xPAM_A c$$$ c$$$ xmi = ( XPP + BETA*YPP - BETA*ALFA )/(1+BETA**2) c$$$ if(xmi.lt.dmin1(xPAM_A,xPAM_B))xmi=dmin1(xPAM_A,xPAM_B) c$$$ if(xmi.gt.dmax1(xPAM_A,xPAM_B))xmi=dmax1(xPAM_A,xPAM_B) c$$$ ymi = ALFA + BETA * xmi c$$$ RE = resyPAM c$$$ c$$$ endif c$$$ c$$$ distance= c$$$ $ ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2 c$$$ distance=dsqrt(distance) c$$$ c$$$c$$$ print*,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b c$$$c$$$ $ ,' --- distance_to --- ',xpp,ypp c$$$c$$$ print*,' resolution ',re c$$$ c$$$ c$$$* ---------------------- c$$$ elseif( c$$$ + xPAM.ne.0.and. c$$$ + yPAM.ne.0.and. c$$$ + zPAM.ne.0.and. c$$$ + xPAM_A.eq.0.and. c$$$ + yPAM_A.eq.0.and. c$$$ + zPAM_A.eq.0.and. c$$$ + xPAM_B.eq.0.and. c$$$ + yPAM_B.eq.0.and. c$$$ + zPAM_B.eq.0.and. c$$$ + .true.)then c$$$* ---------------------- c$$$* DISTANCE TO --- COUPLE c$$$* ---------------------- c$$$ c$$$ distance= c$$$ $ ((xPAM-XPP)/resxPAM)**2 c$$$ $ + c$$$ $ ((yPAM-YPP)/resyPAM)**2 c$$$ distance=dsqrt(distance) c$$$ c$$$c$$$ print*,xPAM,yPAM,zPAM c$$$c$$$ $ ,' --- distance_to --- ',xpp,ypp c$$$c$$$ print*,' resolution ',resxPAM,resyPAM c$$$ c$$$ else c$$$ c$$$ print* c$$$ $ ,' function distance_to ---> wrong usage!!!' c$$$ print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM c$$$ print*,' xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b ' c$$$ $ ,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b c$$$ endif c$$$ c$$$ distance_to = sngl(distance) c$$$ c$$$ return c$$$ end c$$$ c$$$******************************************************************************** c$$$******************************************************************************** c$$$******************************************************************************** c$$$******************************************************************************** c$$$ c$$$ subroutine whichsensor(nplPAM,xPAM,yPAM,ladder,sensor) c$$$* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * c$$$* Given the view and the (xPAM,yPAM) coordinates (in the PAMELA c$$$* reference system), it returns the ladder and the sensor c$$$* which the point belongs to. c$$$* c$$$* The method to assign a point to a sensor consists in c$$$* - calculating the sum of the distances between the point c$$$* and the sensor edges c$$$* - requiring that it is less-equal than (SiDimX+SiDimY) c$$$* (should be strictly equal actually...) c$$$* c$$$* NB -- SiDimX and SiDimY are not the dimentions of the SENSITIVE volume c$$$* but of the whole silicon sensor c$$$* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * c$$$ include '../common/commontracker.f' c$$$ include '../common/common_align.f' c$$$ c$$$ integer ladder,sensor,viewx,viewy c$$$ real c1(4),c2(4),c3(4) c$$$ data c1/1.,0.,0.,1./ c$$$ data c2/1.,-1.,-1.,1./ c$$$ data c3/1.,1.,0.,0./ c$$$ real*8 yvvv,xvvv c$$$ double precision xi,yi,zi c$$$ double precision xrt,yrt,zrt c$$$ real AA,BB c$$$ real yvv(4),xvv(4) c$$$ c$$$* tollerance to consider the track inside the sensitive area c$$$ real ptoll c$$$ data ptoll/150./ !um c$$$ c$$$ external nviewx,nviewy,acoordsi,dcoord c$$$ c$$$ c$$$ c$$$ nplpt = nplPAM !plane c$$$ viewx = nviewx(nplpt) c$$$ viewy = nviewy(nplpt) c$$$ c$$$ do il=1,nladders_view c$$$ do is=1,2 c$$$ c$$$ do iv=1,4 !loop on sensor vertexes c$$$ stripx = (il-c1(iv))*1024 + c1(iv) + c2(iv)*3 c$$$ stripy = (il-c3(iv))*1024 + c3(iv) c$$$c------------------------------------------------------------------------ c$$$c (xi,yi,zi) = mechanical coordinates in the silicon sensor frame c$$$c------------------------------------------------------------------------ c$$$ xi = acoordsi(stripx,viewx) c$$$ yi = acoordsi(stripy,viewy) c$$$ zi = 0. c$$$c------------------------------------------------------------------------ c$$$c (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame c$$$c------------------------------------------------------------------------ c$$$c N.B. I convert angles from microradiants to radiants c$$$ xrt = xi c$$$ $ - omega(nplpt,il,is)*yi c$$$ $ + gamma(nplpt,il,is)*zi c$$$ $ + dx(nplpt,il,is) c$$$ yrt = omega(nplpt,il,is)*xi c$$$ $ + yi c$$$ $ - beta(nplpt,il,is)*zi c$$$ $ + dy(nplpt,il,is) c$$$ zrt = -gamma(nplpt,il,is)*xi c$$$ $ + beta(nplpt,il,is)*yi c$$$ $ + zi c$$$ $ + dz(nplpt,il,is) c$$$c------------------------------------------------------------------------ c$$$c measured coordinates (in cm) in PAMELA reference system c$$$c------------------------------------------------------------------------ c$$$ yvvv = dcoord(yrt,viewy,il,is) / 1.d4 c$$$ xvvv = dcoord(xrt,viewx,il,is) / 1.d4 c$$$ c$$$ yvv(iv)=sngl(yvvv) c$$$ xvv(iv)=sngl(xvvv) c$$$c print*,'LADDER ',il,' SENSOR ',is,' vertexes >> ' c$$$c $ ,iv,xvv(iv),yvv(iv) c$$$ enddo !end loop on sensor vertexes c$$$ c$$$ dtot=0. c$$$ do iside=1,4,2 !loop on sensor edges X c$$$ iv1=iside c$$$ iv2=mod(iside,4)+1 c$$$* straight line passing trhough two consecutive vertexes c$$$ AA = (yvv(iv1)-yvv(iv2))/(xvv(iv1)-xvv(iv2)) c$$$ BB = yvv(iv1) - AA*xvv(iv1) c$$$* point along the straight line closer to the track c$$$ xoo = (xPAM+AA*yPAM-AA*BB)/(1+AA**2) c$$$ yoo = AA*xoo + BB c$$$* sum of the distances c$$$ dtot = dtot + c$$$ $ sqrt((xPAM-xoo)**2+(yPAM-yoo)**2) c$$$ enddo !end loop on sensor edges c$$$ do iside=2,4,2 !loop on sensor edges Y c$$$ iv1=iside c$$$ iv2=mod(iside,4)+1 c$$$* straight line passing trhough two consecutive vertexes c$$$ AA = (xvv(iv1)-xvv(iv2))/(yvv(iv1)-yvv(iv2)) c$$$ BB = xvv(iv1) - AA*yvv(iv1) c$$$* point along the straight line closer to the track c$$$ yoo = (yPAM+AA*xPAM-AA*BB)/(1+AA**2) c$$$ xoo = AA*yoo + BB c$$$* sum of the distances c$$$ dtot = dtot + c$$$ $ sqrt((xPAM-xoo)**2+(yPAM-yoo)**2) c$$$ enddo !end loop on sensor edges c$$$ c$$$ c$$$* half-perimeter of sensitive area c$$$ Perim = c$$$ $ SiDimX - edgeX_l - edgeX_r c$$$ $ +SiDimY - edgeY_l - edgeY_r c$$$ Perim = (Perim + ptoll)/1.e4 c$$$ if(dtot.le.Perim)goto 100 c$$$ c$$$ c$$$ enddo c$$$ enddo c$$$ c$$$ ladder = 0 c$$$ sensor = 0 c$$$ goto 200 c$$$ c$$$ 100 continue c$$$ ladder = il c$$$ sensor = is c$$$ c$$$ c$$$ 200 return c$$$ end c$$$ c$$$ c$$$ c$$$************************************************************************* c$$$ c$$$ subroutine reverse(v,n,temp) !invert the order of the components of v(n) vector c$$$ c$$$ implicit double precision (A-H,O-Z) c$$$ c$$$ dimension v(*) c$$$ dimension temp(*) c$$$ integer i,n c$$$ c$$$ do i=1,n c$$$ temp(i)=v(n+1-i) c$$$ enddo c$$$ c$$$ do i=1,n c$$$ v(i)=temp(i) c$$$ enddo c$$$ c$$$ return c$$$ end c$$$ c$$$************************************************************************* c$$$************************************************************************* c$$$************************************************************************* c$$$************************************************************************* c$$$************************************************************************* c$$$************************************************************************* c$$$************************************************************************* c$$$************************************************************************* c$$$************************************************************************* c$$$************************************************************************* c$$$************************************************************************* c$$$************************************************************************* c$$$************************************************************************* c$$$************************************************************************* c$$$************************************************************************* c$$$ integer function ip_cp(id) c$$$* c$$$* given the couple id, c$$$* it returns the plane number c$$$* c$$$ include '../common/commontracker.f' c$$$c include '../common/common_analysis.f' c$$$ include '../common/common_momanhough.f' c$$$ c$$$ ip_cp=0 c$$$ ncpp=0 c$$$ do ip=1,nplanes c$$$ ncpp=ncpp+ncp_plane(ip) c$$$ if(ncpp.ge.abs(id))then c$$$ ip_cp=ip c$$$ goto 100 c$$$ endif c$$$ enddo c$$$ 100 continue c$$$ return c$$$ end c$$$ c$$$ c$$$ integer function is_cp(id) c$$$* c$$$* given the couple id, c$$$* it returns the sensor number c$$$* c$$$ is_cp=0 c$$$ if(id.lt.0)is_cp=1 c$$$ if(id.gt.0)is_cp=2 c$$$ if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!' c$$$ c$$$ return c$$$ end c$$$ c$$$ c$$$ integer function icp_cp(id) c$$$* c$$$* given the couple id, c$$$* it returns the id number ON THE PLANE c$$$* c$$$ include '../common/commontracker.f' c$$$c include '../common/common_analysis.f' c$$$ include '../common/common_momanhough.f' c$$$ c$$$ icp_cp=0 c$$$ c$$$ ncpp=0 c$$$ do ip=1,nplanes c$$$ ncppold=ncpp c$$$ ncpp=ncpp+ncp_plane(ip) c$$$ if(ncpp.ge.abs(id))then c$$$ icp_cp=abs(id)-ncppold c$$$ goto 100 c$$$ endif c$$$ enddo c$$$ 100 continue c$$$ return c$$$ end c$$$ c$$$ c$$$ c$$$ integer function id_cp(ip,icp,is) c$$$* c$$$* given a plane, a couple and the sensor c$$$* it returns the absolute couple id c$$$* negative if sensor =1 c$$$* positive if sensor =2 c$$$* c$$$ include '../common/commontracker.f' c$$$c include '../common/calib.f' c$$$c include '../common/level1.f' c$$$c include '../common/common_analysis.f' c$$$ include '../common/common_momanhough.f' c$$$ c$$$ id_cp=0 c$$$ c$$$ if(ip.gt.1)then c$$$ do i=1,ip-1 c$$$ id_cp = id_cp + ncp_plane(i) c$$$ enddo c$$$ endif c$$$ c$$$ id_cp = id_cp + icp c$$$ c$$$ if(is.eq.1) id_cp = -id_cp c$$$ c$$$ return c$$$ end c$$$ c$$$ c$$$ c$$$ c$$$************************************************************************* c$$$************************************************************************* c$$$************************************************************************* c$$$************************************************************************* c$$$************************************************************************* c$$$************************************************************************* c$$$ subroutine book_debug c$$$ c$$$ include '../common/commontracker.f' c$$$c include '../common/calib.f' c$$$c include '../common/level1.f' c$$$c include '../common/common_analysis.f' c$$$ include '../common/common_momanhough.f' c$$$ include '../common/common_level2debug.f' c$$$ c$$$c include '../common/common_mini.f' c$$$ c$$$ character*35 block1,block2,block3,block4 c$$$ $ ,block5,block6 c$$$ c$$$ c$$$* * * * * * * * * * * * * * * * * * * * * * * * c$$$* HOUGH TRANSFORM PARAMETERS c$$$c call HBOOK1(1001 c$$$c $ ,'y' c$$$c $ ,3000,-70.,70.,0.) c$$$c call HBOOK1(1002 c$$$c $ ,'tg thyz' c$$$c $ ,300,-1.,1.,0.) c$$$ c$$$ call HBOOK2(1003 c$$$ $ ,'y vs tg thyz' c$$$ $ ,300,-1.,1. !x c$$$ $ ,3000,-70.,70.,0.) !y c$$$ c$$$ call HBOOK1(1004 c$$$ $ ,'Dy' c$$$ $ ,100,0.,2.,0.) c$$$ c$$$ call HBOOK1(1005 c$$$ $ ,'D thyz' c$$$ $ ,100,0.,.05,0.) c$$$ c$$$ c$$$ c$$$* DEBUG ntuple: c$$$ call HBNT(ntp_level2+1,'LEVEL2',' ') c$$$ c$$$ call HBNAME(ntp_level2+1,'EVENT',good2_nt, c$$$ $ 'GOOD2:L,NEV2:I') c$$$ c$$$ 411 format('NDBLT:I::[0,',I5,']') c$$$ write(block1,411) ndblt_max_nt c$$$ call HBNAME(ntp_level2+1,'HOUGH YZ',ndblt_nt, c$$$ $ block1//' c$$$ $ ,ALFAYZ1(NDBLT):R c$$$ $ ,ALFAYZ2(NDBLT):R c$$$ $ ') c$$$ c$$$ 412 format('NTRPT:I::[0,',I5,']') c$$$ write(block2,412) ntrpt_max_nt c$$$ call HBNAME(ntp_level2+1,'HOUGH XZ',NTRPT_nt, c$$$ $ block2//' c$$$ $ ,ALFAXZ1(NTRPT):R c$$$ $ ,ALFAXZ2(NTRPT):R c$$$ $ ,ALFAXZ3(NTRPT):R c$$$ $ ') c$$$ c$$$ c$$$ 413 format('NCLOUDS_YZ:I::[0,',I4,']') c$$$ 414 format('DB_CLOUD(',I4,'):I') c$$$ write(block3,413) ncloyz_max c$$$ write(block4,414) ndblt_max_nt c$$$ call HBNAME(ntp_level2+1,'CLOUD YZ',NCLOUDS_YZ, c$$$ $ block3//' c$$$ $ ,ALFAYZ1_AV(NCLOUDS_YZ):R c$$$ $ ,ALFAYZ2_AV(NCLOUDS_YZ):R c$$$ $ ,PTCLOUD_YZ(NCLOUDS_YZ):I c$$$ $ ,'//block4 c$$$ $ ) c$$$ c$$$ 415 format('NCLOUDS_XZ:I::[0,',I4,']') c$$$ 416 format('TR_CLOUD(',I5,'):I') c$$$ write(block5,415) ncloxz_max c$$$ write(block6,416) ntrpt_max_nt c$$$ call HBNAME(ntp_level2+1,'CLOUD XZ',NCLOUDS_XZ, c$$$ $ block5//' c$$$ $ ,ALFAXZ1_AV(NCLOUDS_XZ):R c$$$ $ ,ALFAXZ2_AV(NCLOUDS_XZ):R c$$$ $ ,ALFAXZ3_AV(NCLOUDS_XZ):R c$$$ $ ,PTCLOUD_XZ(NCLOUDS_XZ):I c$$$ $ ,'//block6 c$$$ $ ) c$$$ c$$$ c$$$ return c$$$ end c$$$***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...*** c$$$* c$$$* c$$$* c$$$* c$$$* c$$$* c$$$* c$$$* c$$$* c$$$***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...*** c$$$ subroutine book_level2 c$$$ c$$$ include '../common/commontracker.f' c$$$ include '../common/common_momanhough.f' c$$$ include '../common/level2.f' c$$$ c$$$ character*35 block1,block2 c$$$ c$$$c print*,'__________ booking LEVEL2 n-tuple __________' c$$$ c$$$* LEVEL1 ntuple: c$$$ call HBNT(ntp_level2,'LEVEL2',' ') c$$$ c$$$ call HBNAME(ntp_level2,'EVENT',good2, c$$$ $ 'GOOD2:L,NEV2:I') c$$$ c$$$# ifndef TEST2003 c$$$ c$$$ call HBNAME(ntp_level2,'CPU',pkt_type c$$$ $ ,'PKT_TYPE:I::[0,50] c$$$ $ ,PKT_NUM:I c$$$ $ ,OBT:I c$$$ $ ,WHICH_CALIB:I::[0,50]') c$$$ c$$$# endif c$$$ c$$$ 417 format('NTRK:I::[0,',I4,']') c$$$ 418 format(',IMAGE(NTRK):I::[0,',I4,']') c$$$ write(block1,417)NTRKMAX c$$$ write(block2,418)NTRKMAX c$$$ call HBNAME(ntp_level2,'TRACKS',NTRK, c$$$ $ block1// c$$$ $ block2//' c$$$ $ ,XM(6,NTRK):R c$$$ $ ,YM(6,NTRK):R c$$$ $ ,ZM(6,NTRK):R c$$$ $ ,RESX(6,NTRK):R c$$$ $ ,RESY(6,NTRK):R c$$$ $ ,AL(5,NTRK):R c$$$ $ ,COVAL(5,5,NTRK):R c$$$ $ ,CHI2(NTRK):R c$$$ $ ,XGOOD(6,NTRK):I::[0,1] c$$$ $ ,YGOOD(6,NTRK):I::[0,1] c$$$ $ ,XV(6,NTRK):R c$$$ $ ,YV(6,NTRK):R c$$$ $ ,ZV(6,NTRK):R c$$$ $ ,AXV(6,NTRK):R c$$$ $ ,AYV(6,NTRK):R c$$$ $ ,DEDXP(6,NTRK):R c$$$ $ ') c$$$ c$$$ c$$$ call HBNAME(ntp_level2,'SINGLETS',nclsx, c$$$ $ 'NCLSX(6):I,NCLSY(6):I') c$$$ c$$$ return c$$$ end c$$$ c$$$* **************************************************** c$$$ c$$$ subroutine init_level2 c$$$ c$$$ include '../common/commontracker.f' c$$$ include '../common/level2.f' c$$$ include '../common/level1.f' c$$$ c$$$ c$$$ good2 = .false. c$$$ nev2 = nev1 c$$$ c$$$# ifndef TEST2003 c$$$ c$$$ pkt_type = pkt_type1 c$$$ pkt_num = pkt_num1 c$$$ obt = obt1 c$$$ which_calib = which_calib1 c$$$ c$$$# endif c$$$ c$$$ NTRK = 0 c$$$ do it=1,NTRACKSMAX c$$$ IMAGE(IT)=0 c$$$ CHI2_NT(IT) = -100000. c$$$ do ip=1,nplanes c$$$ XM_nt(IP,IT) = 0 c$$$ YM_nt(IP,IT) = 0 c$$$ ZM_nt(IP,IT) = 0 c$$$ RESX_nt(IP,IT) = 0 c$$$ RESY_nt(IP,IT) = 0 c$$$ XGOOD_nt(IP,IT) = 0 c$$$ YGOOD_nt(IP,IT) = 0 c$$$ enddo c$$$ do ipa=1,5 c$$$ AL_nt(IPA,IT) = 0 c$$$ do ipaa=1,5 c$$$ coval(ipa,ipaa,IT)=0 c$$$ enddo c$$$ enddo c$$$ enddo c$$$ c$$$ do ip=1,nplanes c$$$ nclsx(ip)=0 c$$$ nclsy(ip)=0 c$$$ enddo c$$$ c$$$ c$$$ end c$$$ c$$$* ------------------------------------------------------- c$$$* This routine fills the ntr-th element of the variables c$$$* inside the level2_tracks common, which correspond c$$$* to the ntr-th track info. c$$$* ------------------------------------------------------- c$$$ c$$$ subroutine fill_level2_tracks(ntr) c$$$ c$$$ include '../common/commontracker.f' c$$$ include '../common/level2.f' c$$$ include '../common/common_mini_2.f' c$$$ c$$$c print*,'TRACK ',ntr c$$$ c$$$ good2=.true. c$$$ chi2_nt(ntr) = sngl(chi2) c$$$c print*,chi2_nt(ntr) c$$$ do i=1,5 c$$$ al_nt(i,ntr) = sngl(al(i)) c$$$ do j=1,5 c$$$ coval(i,j,ntr) = sngl(cov(i,j)) c$$$ enddo c$$$c print*,al_nt(i,ntr) c$$$ enddo c$$$ do ip=1,nplanes ! loop on planes c$$$ xgood_nt(ip,ntr) = int(xgood(ip)) c$$$ ygood_nt(ip,ntr) = int(ygood(ip)) c$$$ xm_nt(ip,ntr) = sngl(xm(ip)) c$$$ ym_nt(ip,ntr) = sngl(ym(ip)) c$$$ zm_nt(ip,ntr) = sngl(zm(ip)) c$$$ RESX_nt(IP,ntr) = sngl(resx(ip)) c$$$ RESY_nt(IP,ntr) = sngl(resy(ip)) c$$$ xv_nt(ip,ntr) = sngl(xv(ip)) c$$$ yv_nt(ip,ntr) = sngl(yv(ip)) c$$$ zv_nt(ip,ntr) = sngl(zv(ip)) c$$$ axv_nt(ip,ntr) = sngl(axv(ip)) c$$$ ayv_nt(ip,ntr) = sngl(ayv(ip)) c$$$ dedxp(ip,ntr) = sngl(dedxtrk(ip)) c$$$c print*,ip,'|',xm_nt(ip,ntr),ym_nt(ip,ntr),zm_nt(ip,ntr) c$$$ enddo c$$$ c$$$ end c$$$ c$$$*************************************************** c$$$* * c$$$* * c$$$* * c$$$* * c$$$* * c$$$* * c$$$************************************************** c$$$ c$$$ subroutine cl_to_couples(iflag) c$$$ c$$$ include '../common/commontracker.f' c$$$ include '../common/common_momanhough.f' c$$$ include '../common/momanhough_init.f' c$$$ include '../common/calib.f' c$$$ include '../common/level1.f' c$$$ c$$$ logical DEBUG c$$$ common/dbg/DEBUG c$$$ c$$$* output flag c$$$* -------------- c$$$* 0 = good event c$$$* 1 = bad event c$$$* -------------- c$$$ integer iflag c$$$ c$$$ integer badseed,badcl c$$$ c$$$* init variables c$$$ ncp_tot=0 c$$$ do ip=1,nplanes c$$$ do ico=1,ncouplemax c$$$ clx(ip,ico)=0 c$$$ cly(ip,ico)=0 c$$$ enddo c$$$ ncp_plane(ip)=0 c$$$ do icl=1,nclstrmax_level2 c$$$ cls(ip,icl)=1 c$$$ enddo c$$$ ncls(ip)=0 c$$$ enddo c$$$ do icl=1,nclstrmax_level2 c$$$ cl_single(icl)=1 c$$$ cl_good(icl)=0 c$$$ enddo c$$$ c$$$* start association c$$$ ncouples=0 c$$$ do icx=1,nclstr1 !loop on cluster (X) c$$$ if(mod(VIEW(icx),2).eq.1)goto 10 c$$$ c$$$* ---------------------------------------------------- c$$$* cut on charge (X VIEW) c$$$ if(dedx(icx).lt.dedx_x_min)then c$$$ cl_single(icx)=0 c$$$ goto 10 c$$$ endif c$$$* cut BAD (X VIEW) c$$$ badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx))) c$$$ ifirst=INDSTART(icx) c$$$ if(icx.ne.nclstr1) then c$$$ ilast=INDSTART(icx+1)-1 c$$$ else c$$$ ilast=TOTCLLENGTH c$$$ endif c$$$ badcl=badseed c$$$ do igood=-ngoodstr,ngoodstr c$$$ ibad=1 c$$$ if((INDMAX(icx)+igood).gt.ifirst.and. c$$$ $ (INDMAX(icx)+igood).lt.ilast.and. c$$$ $ .true.)then c$$$ ibad=BAD(VIEW(icx), c$$$ $ nvk(MAXS(icx)+igood), c$$$ $ nst(MAXS(icx)+igood)) c$$$ endif c$$$ badcl=badcl*ibad c$$$ enddo c$$$c if(badcl.eq.0)then c$$$c cl_single(icx)=0 c$$$c goto 10 c$$$c endif c$$$* ---------------------------------------------------- c$$$ c$$$ cl_good(icx)=1 c$$$ nplx=npl(VIEW(icx)) c$$$ nldx=nld(MAXS(icx),VIEW(icx)) c$$$ c$$$ do icy=1,nclstr1 !loop on cluster (Y) c$$$ if(mod(VIEW(icy),2).eq.0)goto 20 c$$$ c$$$* ---------------------------------------------------- c$$$* cut on charge (Y VIEW) c$$$ if(dedx(icy).lt.dedx_y_min)then c$$$ cl_single(icy)=0 c$$$ goto 20 c$$$ endif c$$$* cut BAD (Y VIEW) c$$$ badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy))) c$$$ ifirst=INDSTART(icy) c$$$ if(icy.ne.nclstr1) then c$$$ ilast=INDSTART(icy+1)-1 c$$$ else c$$$ ilast=TOTCLLENGTH c$$$ endif c$$$ badcl=badseed c$$$ do igood=-ngoodstr,ngoodstr c$$$ ibad=1 c$$$ if((INDMAX(icy)+igood).gt.ifirst.and. c$$$ $ (INDMAX(icy)+igood).lt.ilast.and. c$$$ $ .true.) c$$$ $ ibad=BAD(VIEW(icy), c$$$ $ nvk(MAXS(icy)+igood), c$$$ $ nst(MAXS(icy)+igood)) c$$$ badcl=badcl*ibad c$$$ enddo c$$$c if(badcl.eq.0)then c$$$c cl_single(icy)=0 c$$$c goto 20 c$$$c endif c$$$* ---------------------------------------------------- c$$$ c$$$ c$$$ cl_good(icy)=1 c$$$ nply=npl(VIEW(icy)) c$$$ nldy=nld(MAXS(icy),VIEW(icy)) c$$$ c$$$* ---------------------------------------------- c$$$* CONDITION TO FORM A COUPLE c$$$* ---------------------------------------------- c$$$* geometrical consistency (same plane and ladder) c$$$ if(nply.eq.nplx.and.nldy.eq.nldx)then c$$$* charge correlation c$$$ ddd=(dedx(icy) c$$$ $ -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx)) c$$$ ddd=ddd/sqrt(kch(nplx,nldx)**2+1) c$$$ cut=chcut*sch(nplx,nldx) c$$$ if(abs(ddd).gt.cut)goto 20 !charge not consistent c$$$ c$$$ c$$$* ------------------> COUPLE <------------------ c$$$* check to do not overflow vector dimentions c$$$ if(ncp_plane(nplx).gt.ncouplemax)then c$$$ if(DEBUG)print*, c$$$ $ ' ** warning ** number of identified'// c$$$ $ ' couples on plane ',nplx, c$$$ $ ' exceeds vector dimention'// c$$$ $ ' ( ',ncouplemax,' )' c$$$c good2=.false. c$$$c goto 880 !fill ntp and go to next event c$$$ iflag=1 c$$$ return c$$$ endif c$$$ c$$$ if(ncp_plane(nplx).eq.ncouplemax)then c$$$ if(DEBUG)print*, c$$$ $ '** warning ** number of identified '// c$$$ $ 'couples on plane ',nplx, c$$$ $ 'exceeds vector dimention ' c$$$ $ ,'( ',ncouplemax,' )' c$$$c good2=.false. c$$$c goto 880 !fill ntp and go to next event c$$$ iflag=1 c$$$ return c$$$ endif c$$$ c$$$ ncp_plane(nplx) = ncp_plane(nplx) + 1 c$$$ clx(nplx,ncp_plane(nplx))=icx c$$$ cly(nply,ncp_plane(nplx))=icy c$$$ cl_single(icx)=0 c$$$ cl_single(icy)=0 c$$$ endif c$$$* ---------------------------------------------- c$$$ c$$$ 20 continue c$$$ enddo !end loop on clusters(Y) c$$$ c$$$ 10 continue c$$$ enddo !end loop on clusters(X) c$$$ c$$$ c$$$ do icl=1,nclstr1 c$$$ if(cl_single(icl).eq.1)then c$$$ ip=npl(VIEW(icl)) c$$$ ncls(ip)=ncls(ip)+1 c$$$ cls(ip,ncls(ip))=icl c$$$ endif c$$$ enddo c$$$ c$$$ c$$$ if(DEBUG)then c$$$ print*,'clusters ',nclstr1 c$$$ print*,'good ',(cl_good(i),i=1,nclstr1) c$$$ print*,'singles ',(cl_single(i),i=1,nclstr1) c$$$ print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes) c$$$ endif c$$$ c$$$ do ip=1,6 c$$$ ncp_tot=ncp_tot+ncp_plane(ip) c$$$ enddo c$$$c if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!) c$$$ c$$$ if(ncp_tot.gt.ncp_max)then c$$$ if(DEBUG)print*, c$$$ $ '** warning ** number of identified '// c$$$ $ 'couples exceeds upper limit for Hough tr. ' c$$$ $ ,'( ',ncp_max,' )' c$$$c good2=.false. c$$$c goto 880 !fill ntp and go to next event c$$$ iflag=1 c$$$ return c$$$ endif c$$$ c$$$ return c$$$ end c$$$ c$$$ c$$$*************************************************** c$$$* * c$$$* * c$$$* * c$$$* * c$$$* * c$$$* * c$$$************************************************** c$$$ c$$$ subroutine cp_to_doubtrip(iflag) c$$$ c$$$ include '../common/commontracker.f' c$$$ include '../common/common_momanhough.f' c$$$ include '../common/momanhough_init.f' c$$$ include '../common/common_xyzPAM.f' c$$$ include '../common/calib.f' c$$$ include '../common/level1.f' c$$$ c$$$ logical DEBUG c$$$ common/dbg/DEBUG c$$$ c$$$* output flag c$$$* -------------- c$$$* 0 = good event c$$$* 1 = bad event c$$$* -------------- c$$$ integer iflag c$$$ c$$$ c$$$* ----------------------------- c$$$* DOUBLETS/TRIPLETS coordinates c$$$c double precision xm1,ym1,zm1 c$$$c double precision xm2,ym2,zm2 c$$$c double precision xm3,ym3,zm3 c$$$ c$$$ real xm1,ym1,zm1 c$$$ real xm2,ym2,zm2 c$$$ real xm3,ym3,zm3 c$$$* ----------------------------- c$$$* variable needed for tricircle: c$$$ real xp(3),zp(3)!TRIPLETS coordinates, to find a circle c$$$ EQUIVALENCE (xm1,xp(1)) c$$$ EQUIVALENCE (xm2,xp(2)) c$$$ EQUIVALENCE (xm3,xp(3)) c$$$ EQUIVALENCE (zm1,zp(1)) c$$$ EQUIVALENCE (zm2,zp(2)) c$$$ EQUIVALENCE (zm3,zp(3)) c$$$ real angp(3),resp(3),chi c$$$ real xc,zc,radius c$$$* ----------------------------- c$$$ c$$$ c$$$ c$$$ ndblt=0 !number of doublets c$$$ ntrpt=0 !number of triplets c$$$ c$$$ do ip1=1,(nplanes-1) !loop on planes - COPPIA 1 c$$$ do is1=1,2 !loop on sensors - COPPIA 1 c$$$ c$$$ do icp1=1,ncp_plane(ip1) !loop on COPPIA 1 c$$$ icx1=clx(ip1,icp1) c$$$ icy1=cly(ip1,icp1) c$$$ call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.) c$$$ xm1=xPAM c$$$ ym1=yPAM c$$$ zm1=zPAM c$$$c print*,'***',is1,xm1,ym1,zm1 c$$$ do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2 c$$$ do is2=1,2 !loop on sensors -ndblt COPPIA 2 c$$$ c$$$ do icp2=1,ncp_plane(ip2) !loop on COPPIA 2 c$$$ icx2=clx(ip2,icp2) c$$$ icy2=cly(ip2,icp2) c$$$ call xyz_PAM c$$$ $ (icx2,icy2,is2,'COG2','COG2',0.,0.) c$$$ xm2=xPAM c$$$ ym2=yPAM c$$$ zm2=zPAM c$$$ c$$$* - - - - - - - - - - - - - - - - - - - - - - - - - - - - c$$$* track parameters on Y VIEW c$$$* (2 couples needed) c$$$* - - - - - - - - - - - - - - - - - - - - - - - - - - - - c$$$ if(ndblt.eq.ndblt_max)then c$$$ if(DEBUG)print*, c$$$ $ '** warning ** number of identified '// c$$$ $ 'doublets exceeds vector dimention ' c$$$ $ ,'( ',ndblt_max,' )' c$$$c good2=.false. c$$$c goto 880 !fill ntp and go to next event c$$$ iflag=1 c$$$ return c$$$ endif c$$$ ndblt = ndblt + 1 c$$$* store doublet info c$$$ cpyz1(ndblt)=id_cp(ip1,icp1,is1) c$$$ cpyz2(ndblt)=id_cp(ip2,icp2,is2) c$$$* tg(th_yz) c$$$ alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2) c$$$* y0 (cm) c$$$ alfayz1(ndblt)=alfayz2(ndblt)*(zini-zm1)+ym1 c$$$ c$$$**** -----------------------------------------------**** c$$$**** reject non phisical couples **** c$$$**** -----------------------------------------------**** c$$$ if( c$$$ $ abs(alfayz2(ndblt)).gt.alfyz2_max c$$$ $ .or. c$$$ $ abs(alfayz1(ndblt)).gt.alfyz1_max c$$$ $ )ndblt = ndblt-1 c$$$ c$$$c$$$ if(iev.eq.33)then c$$$c$$$ print*,'********* ',ndblt,' -- ',icp1,icp2,is1,is2 c$$$c$$$ $ ,' || ',icx1,icy1,icx2,icy2 c$$$c$$$ $ ,' || ',xm1,ym1,xm2,ym2 c$$$c$$$ $ ,' || ',alfayz2(ndblt),alfayz1(ndblt) c$$$c$$$ endif c$$$c$$$ c$$$* - - - - - - - - - - - - - - - - - - - - - - - - - - - - c$$$* track parameters on Y VIEW - end c$$$* - - - - - - - - - - - - - - - - - - - - - - - - - - - - c$$$ c$$$ c$$$ if(ip2.eq.nplanes)goto 30 !no possible combination with 3 couples c$$$ do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3 c$$$ do is3=1,2 !loop on sensors - COPPIA 3 c$$$ c$$$ do icp3=1,ncp_plane(ip3) !loop on COPPIA 3 c$$$ icx3=clx(ip3,icp3) c$$$ icy3=cly(ip3,icp3) c$$$ call xyz_PAM c$$$ $ (icx3,icy3,is3,'COG2','COG2',0.,0.) c$$$ xm3=xPAM c$$$ ym3=yPAM c$$$ zm3=zPAM c$$$* find the circle passing through the three points c$$$ call tricircle(3,xp,zp,angp,resp,chi c$$$ $ ,xc,zc,radius,iflag) c$$$c print*,xc,zc,radius c$$$* the circle must intersect the reference plane c$$$ if( c$$$c $ (xc.le.-1.*xclimit.or. c$$$c $ xc.ge.xclimit).and. c$$$ $ radius**2.ge.(ZINI-zc)**2.and. c$$$ $ iflag.eq.0.and. c$$$ $ .true.)then c$$$ c$$$* - - - - - - - - - - - - - - - - - - - - - - - - - - - - c$$$* track parameters on X VIEW c$$$* (3 couples needed) c$$$* - - - - - - - - - - - - - - - - - - - - - - - - - - - - c$$$ if(ntrpt.eq.ntrpt_max)then c$$$ if(DEBUG)print*, c$$$ $ '** warning ** number of identified '// c$$$ $ 'triplets exceeds vector dimention ' c$$$ $ ,'( ',ntrpt_max,' )' c$$$c good2=.false. c$$$c goto 880 !fill ntp and go to next event c$$$ iflag=1 c$$$ return c$$$ endif c$$$ ntrpt = ntrpt +1 c$$$* store triplet info c$$$ cpxz1(ntrpt)=id_cp(ip1,icp1,is1) c$$$ cpxz2(ntrpt)=id_cp(ip2,icp2,is2) c$$$ cpxz3(ntrpt)=id_cp(ip3,icp3,is3) c$$$ c$$$ if(xc.lt.0)then c$$$*************POSITIVE DEFLECTION c$$$ alfaxz1(ntrpt) = xc+sqrt(radius**2-(ZINI-zc)**2) c$$$ alfaxz2(ntrpt) = (ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2) c$$$ alfaxz3(ntrpt) = 1/radius c$$$ else c$$$*************NEGATIVE DEFLECTION c$$$ alfaxz1(ntrpt) = xc-sqrt(radius**2-(ZINI-zc)**2) c$$$ alfaxz2(ntrpt) = -(ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2) c$$$ alfaxz3(ntrpt) = -1/radius c$$$ endif c$$$ c$$$**** -----------------------------------------------**** c$$$**** reject non phisical triplets **** c$$$**** -----------------------------------------------**** c$$$ if( c$$$ $ abs(alfaxz2(ntrpt)).gt.alfxz2_max c$$$ $ .or. c$$$ $ abs(alfaxz1(ntrpt)).gt.alfxz1_max c$$$ $ )ntrpt = ntrpt-1 c$$$ c$$$ c$$$c print*,alfaxz1(ntrpt),alfaxz2(ntrpt),alfaxz3(ntrpt) c$$$* - - - - - - - - - - - - - - - - - - - - - - - - - - - - c$$$* track parameters on X VIEW - end c$$$* - - - - - - - - - - - - - - - - - - - - - - - - - - - - c$$$ endif c$$$ enddo !end loop on COPPIA 3 c$$$ enddo !end loop on sensors - COPPIA 3 c$$$ enddo !end loop on planes - COPPIA 3 c$$$ 30 continue c$$$ c$$$ 1 enddo !end loop on COPPIA 2 c$$$ enddo !end loop on sensors - COPPIA 2 c$$$ enddo !end loop on planes - COPPIA 2 c$$$ c$$$ enddo !end loop on COPPIA1 c$$$ enddo !end loop on sensors - COPPIA 1 c$$$ enddo !end loop on planes - COPPIA 1 c$$$ c$$$ if(DEBUG)then c$$$ print*,'--- doublets ',ndblt c$$$ print*,'--- triplets ',ntrpt c$$$ endif c$$$ c$$$c goto 880 !ntp fill c$$$ c$$$ c$$$ return c$$$ end c$$$ c$$$ c$$$ c$$$*************************************************** c$$$* * c$$$* * c$$$* * c$$$* * c$$$* * c$$$* * c$$$************************************************** c$$$ c$$$ subroutine doub_to_YZcloud(iflag) c$$$ c$$$ include '../common/commontracker.f' c$$$ include '../common/common_momanhough.f' c$$$ include '../common/momanhough_init.f' c$$$ c$$$ logical DEBUG c$$$ common/dbg/DEBUG c$$$ c$$$* output flag c$$$* -------------- c$$$* 0 = good event c$$$* 1 = bad event c$$$* -------------- c$$$ integer iflag c$$$ c$$$ integer db_used(ndblt_max) c$$$ integer db_temp(ndblt_max) c$$$ integer db_all(ndblt_max) !stores db ID in each cloud c$$$ c$$$ integer hit_plane(nplanes) c$$$ c$$$* mask for used couples c$$$ integer cp_useds1(ncouplemaxtot) ! sensor 1 c$$$ integer cp_useds2(ncouplemaxtot) ! sensor 2 c$$$ c$$$* ****************************************** c$$$* PLANE TO EXCLUDE (for efficiency studies) c$$$* --> provided as input parameter c$$$* (set according to "tracking" notation): c$$$* 1-6 from TOP to BOTTOM c$$$* ****************************************** c$$$ integer ipexclude c$$$ common/plane_exclude/ipexclude c$$$ c$$$ c$$$*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c$$$* classification of DOUBLETS c$$$* according to distance in parameter space c$$$* (cloud = group of points (doublets) in parameter space) c$$$*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c$$$ do idb=1,ndblt c$$$ db_used(idb)=0 c$$$ enddo c$$$ c$$$ distance=0 c$$$ nclouds_yz=0 !number of clouds c$$$ npt_tot=0 c$$$ do idb1=1,ndblt !loop (1) on DOUBLETS c$$$ if(db_used(idb1).eq.1)goto 2228 !db already included in a cloud c$$$* *************************************** c$$$* Checks if a plane should be excluded c$$$* (for efficiency study) c$$$* *************************************** c$$$ ip1=ip_cp(cpyz1(idb1)) c$$$ ip2=ip_cp(cpyz2(idb1)) c$$$ if( c$$$ $ ip1.eq.(nplanes-ipexclude+1).or. c$$$ $ ip2.eq.(nplanes-ipexclude+1).or. c$$$ $ .false.)goto 2228 !exclude DOUBLET from cloud c$$$* *************************************** c$$$ c$$$c print*,'--------------' c$$$c print*,'** ',idb1,' **' c$$$ c$$$ do icp=1,ncp_tot c$$$ cp_useds1(icp)=0 !init c$$$ cp_useds2(icp)=0 !init c$$$ enddo c$$$ do idb=1,ndblt c$$$ db_all(idb)=0 c$$$ enddo c$$$ if(cpyz1(idb1).gt.0)cp_useds2(cpyz1(idb1))=1 c$$$ if(cpyz1(idb1).lt.0)cp_useds1(-cpyz1(idb1))=1 c$$$ if(cpyz2(idb1).gt.0)cp_useds2(cpyz2(idb1))=1 c$$$ if(cpyz2(idb1).lt.0)cp_useds1(-cpyz2(idb1))=1 c$$$ temp1 = alfayz1(idb1) c$$$ temp2 = alfayz2(idb1) c$$$ npt=1 !counter of points in the cloud c$$$ c$$$ db_all(npt) = idb1 c$$$ c$$$ nptloop=1 c$$$ db_temp(1)=idb1 c$$$ c$$$ 88 continue c$$$ c$$$ npv=0 !# new points inlcuded c$$$ do iloop=1,nptloop c$$$ idbref=db_temp(iloop) !local point of reference c$$$ccccc if(db_used(idbref).eq.1)goto 1188 !next c$$$ c$$$ do idb2=1,ndblt !loop (2) on DOUBLETS c$$$ if(idb2.eq.idbref)goto 1118 !next doublet c$$$ if(db_used(idb2).eq.1)goto 1118 c$$$* *************************************** c$$$* Checks if a plane should be excluded c$$$* (for efficiency study) c$$$* *************************************** c$$$ ip1=ip_cp(cpyz1(idb2)) c$$$ ip2=ip_cp(cpyz2(idb2)) c$$$ if( c$$$ $ ip1.eq.(nplanes-ipexclude+1).or. c$$$ $ ip2.eq.(nplanes-ipexclude+1).or. c$$$ $ .false.)goto 1118 !exclude DOUBLET from cloud c$$$* *************************************** c$$$ c$$$ c$$$* doublet distance in parameter space c$$$ distance= c$$$ $ ((alfayz1(idbref)-alfayz1(idb2))/Dalfayz1)**2 c$$$ $ +((alfayz2(idbref)-alfayz2(idb2))/Dalfayz2)**2 c$$$ distance = sqrt(distance) c$$$ c$$$c$$$ if(iev.eq.33)then c$$$c$$$ if(distance.lt.100) c$$$c$$$ $ print*,'********* ',idb1,idbref,idb2,distance c$$$c$$$ if(distance.lt.100) c$$$c$$$ $ print*,'********* ',alfayz1(idbref),alfayz1(idb2) c$$$c$$$ $ ,alfayz2(idbref),alfayz2(idb2) c$$$c$$$ endif c$$$ if(distance.lt.cutdistyz)then c$$$ c$$$c print*,idb1,idb2,distance,' cloud ',nclouds_yz c$$$ if(cpyz1(idb2).gt.0)cp_useds2(cpyz1(idb2))=1 c$$$ if(cpyz1(idb2).lt.0)cp_useds1(-cpyz1(idb2))=1 c$$$ if(cpyz2(idb2).gt.0)cp_useds2(cpyz2(idb2))=1 c$$$ if(cpyz2(idb2).lt.0)cp_useds1(-cpyz2(idb2))=1 c$$$ npt = npt + 1 !counter of points in the cloud c$$$ c$$$ npv = npv +1 c$$$ db_temp(npv) = idb2 c$$$ db_used(idbref) = 1 c$$$ db_used(idb2) = 1 c$$$ c$$$ db_all(npt) = idb2 c$$$ c$$$ temp1 = temp1 + alfayz1(idb2) c$$$ temp2 = temp2 + alfayz2(idb2) c$$$c print*,'* idbref,idb2 ',idbref,idb2 c$$$ endif c$$$ c$$$ 1118 continue c$$$ enddo !end loop (2) on DOUBLETS c$$$ c$$$ 1188 continue c$$$ enddo !end loop on... bo? c$$$ c$$$ nptloop=npv c$$$ if(nptloop.ne.0)goto 88 c$$$ c$$$* ------------------------------------------ c$$$* stores the cloud only if c$$$* 1) it includes a minimum number of REAL couples c$$$* 1bis) it inlcudes a minimum number of doublets c$$$* 2) it is not already stored c$$$* ------------------------------------------ c$$$ do ip=1,nplanes c$$$ hit_plane(ip)=0 c$$$ enddo c$$$ ncpused=0 c$$$ do icp=1,ncp_tot c$$$ if(cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0)then c$$$ ncpused=ncpused+1 c$$$ ip=ip_cp(icp) c$$$ hit_plane(ip)=1 c$$$ endif c$$$ enddo c$$$ nplused=0 c$$$ do ip=1,nplanes c$$$ nplused=nplused+ hit_plane(ip) c$$$ enddo c$$$c print*,'>>>> ',ncpused,npt,nplused c$$$ if(ncpused.lt.ncpyz_min)goto 2228 !next doublet c$$$ if(npt.lt.nptyz_min)goto 2228 !next doublet c$$$ if(nplused.lt.nplyz_min)goto 2228 !next doublet c$$$ c$$$* ~~~~~~~~~~~~~~~~~ c$$$* >>> NEW CLOUD <<< c$$$ c$$$ if(nclouds_yz.ge.ncloyz_max)then c$$$ if(DEBUG)print*, c$$$ $ '** warning ** number of identified '// c$$$ $ 'YZ clouds exceeds vector dimention ' c$$$ $ ,'( ',ncloyz_max,' )' c$$$c good2=.false. c$$$c goto 880 !fill ntp and go to next event c$$$ iflag=1 c$$$ return c$$$ endif c$$$ c$$$ nclouds_yz = nclouds_yz + 1 !increase counter c$$$ alfayz1_av(nclouds_yz) = temp1/npt !store average parameter c$$$ alfayz2_av(nclouds_yz) = temp2/npt ! " c$$$ do icp=1,ncp_tot c$$$ cpcloud_yz(nclouds_yz,icp)= c$$$ $ cp_useds1(icp)+2*cp_useds2(icp) !store cp info c$$$ enddo c$$$ ptcloud_yz(nclouds_yz)=npt c$$$c ptcloud_yz_nt(nclouds_yz)=npt c$$$ do ipt=1,npt c$$$ db_cloud(npt_tot+ipt) = db_all(ipt) c$$$c print*,'>> ',ipt,db_all(ipt) c$$$ enddo c$$$ npt_tot=npt_tot+npt c$$$ if(DEBUG)then c$$$ print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~' c$$$ print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points' c$$$ print*,'- alfayz1 ',alfayz1_av(nclouds_yz) c$$$ print*,'- alfayz2 ',alfayz2_av(nclouds_yz) c$$$ print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot) c$$$ print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot) c$$$ print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes) c$$$ endif c$$$* >>> NEW CLOUD <<< c$$$* ~~~~~~~~~~~~~~~~~ c$$$ 2228 continue c$$$ enddo !end loop (1) on DOUBLETS c$$$ c$$$ c$$$ if(DEBUG)then c$$$ print*,'---------------------- ' c$$$ print*,'Y-Z total clouds ',nclouds_yz c$$$ print*,' ' c$$$ endif c$$$ c$$$ c$$$ return c$$$ end c$$$ c$$$ c$$$ c$$$ c$$$ c$$$*************************************************** c$$$* * c$$$* * c$$$* * c$$$* * c$$$* * c$$$* * c$$$************************************************** c$$$ c$$$ subroutine trip_to_XZcloud(iflag) c$$$ c$$$ include '../common/commontracker.f' c$$$ include '../common/common_momanhough.f' c$$$ include '../common/momanhough_init.f' c$$$ c$$$ logical DEBUG c$$$ common/dbg/DEBUG c$$$ c$$$* output flag c$$$* -------------- c$$$* 0 = good event c$$$* 1 = bad event c$$$* -------------- c$$$ integer iflag c$$$ c$$$ integer tr_used(ntrpt_max) c$$$ integer tr_temp(ntrpt_max) c$$$ integer tr_incl(ntrpt_max) c$$$ integer tr_all(ntrpt_max) !stores tr ID in each cloud c$$$ c$$$ integer hit_plane(nplanes) c$$$ c$$$* mask for used couples c$$$ integer cp_useds1(ncouplemaxtot) ! sensor 1 c$$$ integer cp_useds2(ncouplemaxtot) ! sensor 2 c$$$ c$$$* ****************************************** c$$$* PLANE TO EXCLUDE (for efficiency studies) c$$$* --> provided as input parameter c$$$* (set according to "tracking" notation): c$$$* 1-6 from TOP to BOTTOM c$$$* ****************************************** c$$$ integer ipexclude c$$$ common/plane_exclude/ipexclude c$$$ c$$$ c$$$*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c$$$* classification of TRIPLETS c$$$* according to distance in parameter space c$$$*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c$$$ do itr=1,ntrpt c$$$ tr_used(itr)=0 c$$$ enddo c$$$ c$$$ distance=0 c$$$ nclouds_xz=0 !number of clouds c$$$ npt_tot=0 !total number of selected triplets c$$$ do itr1=1,ntrpt !loop (1) on TRIPLETS c$$$ if(tr_used(itr1).eq.1)goto 22288 !already included in a cloud c$$$* *************************************** c$$$* Checks if a plane should be excluded c$$$* (for efficiency study) c$$$* *************************************** c$$$ ip1=ip_cp(cpxz1(itr1)) c$$$ ip2=ip_cp(cpxz2(itr1)) c$$$ ip3=ip_cp(cpxz3(itr1)) c$$$ if( c$$$ $ ip1.eq.(nplanes-ipexclude+1).or. c$$$ $ ip2.eq.(nplanes-ipexclude+1).or. c$$$ $ ip3.eq.(nplanes-ipexclude+1).or. c$$$ $ .false.)goto 22288 !exclude TRIPLET from cloud c$$$* *************************************** c$$$c print*,'--------------' c$$$c print*,'** ',itr1,' **' c$$$ c$$$ do icp=1,ncp_tot c$$$ cp_useds1(icp)=0 c$$$ cp_useds2(icp)=0 c$$$ enddo c$$$ do itr=1,ntrpt c$$$ tr_all(itr)=0 !list of included triplets c$$$ enddo c$$$ if(cpxz1(itr1).gt.0)cp_useds2(cpxz1(itr1))=1 c$$$ if(cpxz1(itr1).lt.0)cp_useds1(-cpxz1(itr1))=1 c$$$ if(cpxz2(itr1).gt.0)cp_useds2(cpxz2(itr1))=1 c$$$ if(cpxz2(itr1).lt.0)cp_useds1(-cpxz2(itr1))=1 c$$$ if(cpxz3(itr1).gt.0)cp_useds2(cpxz3(itr1))=1 c$$$ if(cpxz3(itr1).lt.0)cp_useds1(-cpxz3(itr1))=1 c$$$ temp1 = alfaxz1(itr1) c$$$ temp2 = alfaxz2(itr1) c$$$ temp3 = alfaxz3(itr1) c$$$ npt=1 !counter of points in the cloud c$$$ c$$$ tr_all(npt) = itr1 c$$$ c$$$ nptloop=1 c$$$c tr_temp(1)=itr1 c$$$ tr_incl(1)=itr1 c$$$ c$$$ 8881 continue c$$$ c$$$c print*,'start search pv ',nptloop c$$$ npv=0 !# new points inlcuded c$$$ do iloop=1,nptloop c$$$c itrref=tr_temp(iloop) !local point of reference c$$$ itrref=tr_incl(iloop) !local point of reference c$$$c print*,'ref ',itrref c$$$ccccc if(tr_used(itrref).eq.1)goto 11888 !next c$$$ do itr2=1,ntrpt !loop (2) on TRIPLETS c$$$ if(itr2.eq.itr1)goto 11188 !next triplet c$$$ if(tr_used(itr2).eq.1)goto 11188 !next triplet c$$$* *************************************** c$$$* Checks if a plane should be excluded c$$$* (for efficiency study) c$$$* *************************************** c$$$ ip1=ip_cp(cpxz1(itr2)) c$$$ ip2=ip_cp(cpxz2(itr2)) c$$$ ip3=ip_cp(cpxz3(itr2)) c$$$ if( c$$$ $ ip1.eq.(nplanes-ipexclude+1).or. c$$$ $ ip2.eq.(nplanes-ipexclude+1).or. c$$$ $ ip3.eq.(nplanes-ipexclude+1).or. c$$$ $ .false.)goto 11188 !exclude TRIPLET from cloud c$$$* *************************************** c$$$ c$$$* triplet distance in parameter space c$$$* solo i due parametri spaziali per il momemnto c$$$ distance= c$$$ $ ((alfaxz1(itrref)-alfaxz1(itr2))/Dalfaxz1)**2 c$$$ $ +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2 c$$$ distance = sqrt(distance) c$$$ c$$$ if(distance.lt.cutdistxz)then c$$$c print*,idb1,idb2,distance,' cloud ',nclouds_yz c$$$ if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1 c$$$ if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1 c$$$ if(cpxz2(itr2).gt.0)cp_useds2(cpxz2(itr2))=1 c$$$ if(cpxz2(itr2).lt.0)cp_useds1(-cpxz2(itr2))=1 c$$$ if(cpxz3(itr2).gt.0)cp_useds2(cpxz3(itr2))=1 c$$$ if(cpxz3(itr2).lt.0)cp_useds1(-cpxz3(itr2))=1 c$$$ npt = npt + 1 !counter of points in the cloud c$$$ c$$$ npv = npv +1 c$$$ tr_temp(npv) = itr2 c$$$ tr_used(itrref) = 1 c$$$ tr_used(itr2) = 1 c$$$ c$$$ tr_all(npt) = itr2 c$$$ c$$$ temp1 = temp1 + alfaxz1(itr2) c$$$ temp2 = temp2 + alfaxz2(itr2) c$$$ temp3 = temp3 + alfaxz3(itr2) c$$$c print*,'* itrref,itr2 ',itrref,itr2,distance c$$$ endif c$$$ c$$$11188 continue c$$$ enddo !end loop (2) on TRIPLETS c$$$ c$$$11888 continue c$$$ enddo !end loop on... bo? c$$$ c$$$ nptloop=npv c$$$ do i=1,npv c$$$ tr_incl(i)=tr_temp(i) c$$$ enddo c$$$ if(nptloop.ne.0)goto 8881 c$$$ c$$$* ------------------------------------------ c$$$* stores the cloud only if c$$$* 1) it includes a minimum number of REAL couples c$$$* 1bis) c$$$* 2) it is not already stored c$$$* ------------------------------------------ c$$$c print*,'check cp_used' c$$$ do ip=1,nplanes c$$$ hit_plane(ip)=0 c$$$ enddo c$$$ ncpused=0 c$$$ do icp=1,ncp_tot c$$$ if(cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0)then c$$$ ncpused=ncpused+1 c$$$ ip=ip_cp(icp) c$$$ hit_plane(ip)=1 c$$$ endif c$$$ enddo c$$$ nplused=0 c$$$ do ip=1,nplanes c$$$ nplused=nplused+ hit_plane(ip) c$$$ enddo c$$$ if(ncpused.lt.ncpxz_min)goto 22288 !next triplet c$$$ if(npt.lt.nptxz_min)goto 22288 !next triplet c$$$ if(nplused.lt.nplxz_min)goto 22288 !next doublet c$$$ c$$$* ~~~~~~~~~~~~~~~~~ c$$$* >>> NEW CLOUD <<< c$$$ if(nclouds_xz.ge.ncloxz_max)then c$$$ if(DEBUG)print*, c$$$ $ '** warning ** number of identified '// c$$$ $ 'XZ clouds exceeds vector dimention ' c$$$ $ ,'( ',ncloxz_max,' )' c$$$c good2=.false. c$$$c goto 880 !fill ntp and go to next event c$$$ iflag=1 c$$$ return c$$$ endif c$$$ nclouds_xz = nclouds_xz + 1 !increase counter c$$$ alfaxz1_av(nclouds_xz) = temp1/npt !store average parameter c$$$ alfaxz2_av(nclouds_xz) = temp2/npt ! " c$$$ alfaxz3_av(nclouds_xz) = temp3/npt ! " c$$$ do icp=1,ncp_tot c$$$ cpcloud_xz(nclouds_xz,icp)= c$$$ $ cp_useds1(icp)+2*cp_useds2(icp) !store cp info c$$$ enddo c$$$ ptcloud_xz(nclouds_xz)=npt c$$$c ptcloud_xz_nt(nclouds_xz)=npt c$$$ do ipt=1,npt c$$$ tr_cloud(npt_tot+ipt) = tr_all(ipt) c$$$c print *,nclouds_xz,ipt,tr_all(ipt),npt_tot+ipt, c$$$c $ tr_cloud(npt_tot+ipt) c$$$ enddo c$$$ npt_tot=npt_tot+npt c$$$ c$$$ if(DEBUG)then c$$$ print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~' c$$$ print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points' c$$$ print*,'- alfaxz1 ',alfaxz1_av(nclouds_xz) c$$$ print*,'- alfaxz2 ',alfaxz2_av(nclouds_xz) c$$$ print*,'- alfaxz3 ',alfaxz3_av(nclouds_xz) c$$$ print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot) c$$$ print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot) c$$$ print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes) c$$$ endif c$$$* >>> NEW CLOUD <<< c$$$* ~~~~~~~~~~~~~~~~~ c$$$22288 continue c$$$ enddo !end loop (1) on DOUBLETS c$$$ c$$$ if(DEBUG)then c$$$ print*,'---------------------- ' c$$$ print*,'X-Z total clouds ',nclouds_xz c$$$ print*,' ' c$$$ endif c$$$ c$$$ c$$$ return c$$$ end c$$$ c$$$ c$$$*************************************************** c$$$* * c$$$* * c$$$* * c$$$* * c$$$* * c$$$* * c$$$************************************************** c$$$ c$$$ subroutine clouds_to_ctrack(iflag) c$$$ c$$$ include '../common/commontracker.f' c$$$ include '../common/common_momanhough.f' c$$$ include '../common/common_xyzPAM.f' c$$$ include '../common/common_mini_2.f' c$$$ include '../common/common_mech.f' c$$$ include '../common/momanhough_init.f' c$$$ c$$$ logical DEBUG c$$$ common/dbg/DEBUG c$$$ c$$$* output flag c$$$* -------------- c$$$* 0 = good event c$$$* 1 = bad event c$$$* -------------- c$$$ integer iflag c$$$ c$$$* ----------------------------------------------------------- c$$$* mask to store (locally) the couples included c$$$* in the intersection bewteen a XZ and YZ cloud c$$$ integer cpintersec(ncouplemaxtot) c$$$* ----------------------------------------------------------- c$$$* list of matching couples in the combination c$$$* between a XZ and YZ cloud c$$$ integer cp_match(nplanes,ncouplemax) c$$$ integer ncp_match(nplanes) c$$$* ----------------------------------------------------------- c$$$ integer hit_plane(nplanes) c$$$* ----------------------------------------------------------- c$$$* variables for track fitting c$$$ double precision AL_INI(5) c$$$ double precision tath c$$$* ----------------------------------------------------------- c$$$c real fitz(nplanes) !z coordinates of the planes in cm c$$$ c$$$ c$$$ ntracks=0 !counter of track candidates c$$$ c$$$ do iyz=1,nclouds_yz !loop on YZ couds c$$$ do ixz=1,nclouds_xz !loop on XZ couds c$$$ c$$$* -------------------------------------------------- c$$$* check of consistency of the clouds c$$$* ---> required a minimum number of matching couples c$$$* the track fit will be performed on the INTERSECTION c$$$* of the two clouds c$$$* -------------------------------------------------- c$$$ do ip=1,nplanes c$$$ hit_plane(ip)=0 c$$$ ncp_match(ip)=0 c$$$ do icpp=1,ncouplemax c$$$ cp_match(ip,icpp)=0 !init couple list c$$$ enddo c$$$ enddo c$$$ ncp_ok=0 c$$$ do icp=1,ncp_tot !loop on couples c$$$* get info on c$$$ cpintersec(icp)=min( c$$$ $ cpcloud_yz(iyz,icp), c$$$ $ cpcloud_xz(ixz,icp)) c$$$ if( c$$$ $ (cpcloud_yz(iyz,icp).eq.1.and.cpcloud_xz(ixz,icp).eq.2).or. c$$$ $ (cpcloud_yz(iyz,icp).eq.2.and.cpcloud_xz(ixz,icp).eq.1).or. c$$$ $ .false.)cpintersec(icp)=0 c$$$ if(cpintersec(icp).ne.0)then c$$$ ncp_ok=ncp_ok+1 c$$$ c$$$ ip=ip_cp(icp) c$$$ hit_plane(ip)=1 c$$$ if(cpintersec(icp).eq.1)then c$$$* 1) only the couple image in sensor 1 matches c$$$ id=-icp c$$$ ncp_match(ip)=ncp_match(ip)+1 c$$$ cp_match(ip,ncp_match(ip))=id c$$$ elseif(cpintersec(icp).eq.2)then c$$$* 2) only the couple image in sensor 2 matches c$$$ id=icp c$$$ ncp_match(ip)=ncp_match(ip)+1 c$$$ cp_match(ip,ncp_match(ip))=id c$$$ else c$$$* 3) both couple images match c$$$ id=icp c$$$ do is=1,2 c$$$ id=-id c$$$ ncp_match(ip)=ncp_match(ip)+1 c$$$ cp_match(ip,ncp_match(ip))=id c$$$ enddo c$$$ endif c$$$ endif !end matching condition c$$$ enddo !end loop on couples c$$$ c$$$ nplused=0 c$$$ do ip=1,nplanes c$$$ nplused=nplused+ hit_plane(ip) c$$$ enddo c$$$ c$$$ if(nplused.lt.nplxz_min)goto 888 !next doublet c$$$ if(ncp_ok.lt.ncpok)goto 888 !next cloud c$$$ c$$$ if(DEBUG)then c$$$ print*,'Combination ',iyz,ixz c$$$ $ ,' db ',ptcloud_yz(iyz) c$$$ $ ,' tr ',ptcloud_xz(ixz) c$$$ $ ,' -----> # matching couples ',ncp_ok c$$$ endif c$$$c$$$ print*,'~~~~~~~~~~~~~~~~~~~~~~~~~' c$$$c$$$ print*,'Configurazione cluster XZ' c$$$c$$$ print*,'1 -- ',(clx(1,i),i=1,ncp_plane(1)) c$$$c$$$ print*,'2 -- ',(clx(2,i),i=1,ncp_plane(1)) c$$$c$$$ print*,'3 -- ',(clx(3,i),i=1,ncp_plane(1)) c$$$c$$$ print*,'4 -- ',(clx(4,i),i=1,ncp_plane(1)) c$$$c$$$ print*,'5 -- ',(clx(5,i),i=1,ncp_plane(1)) c$$$c$$$ print*,'6 -- ',(clx(6,i),i=1,ncp_plane(1)) c$$$c$$$ print*,'Configurazione cluster YZ' c$$$c$$$ print*,'1 -- ',(cly(1,i),i=1,ncp_plane(1)) c$$$c$$$ print*,'2 -- ',(cly(2,i),i=1,ncp_plane(1)) c$$$c$$$ print*,'3 -- ',(cly(3,i),i=1,ncp_plane(1)) c$$$c$$$ print*,'4 -- ',(cly(4,i),i=1,ncp_plane(1)) c$$$c$$$ print*,'5 -- ',(cly(5,i),i=1,ncp_plane(1)) c$$$c$$$ print*,'6 -- ',(cly(6,i),i=1,ncp_plane(1)) c$$$c$$$ print*,'~~~~~~~~~~~~~~~~~~~~~~~~~' c$$$ c$$$* -------> INITIAL GUESS <------- c$$$ AL_INI(1)=dreal(alfaxz1_av(ixz)) c$$$ AL_INI(2)=dreal(alfayz1_av(iyz)) c$$$ AL_INI(4)=datan(dreal(alfayz2_av(iyz)) c$$$ $ /dreal(alfaxz2_av(ixz))) c$$$ tath=-dreal(alfaxz2_av(ixz))/dcos(AL_INI(4)) c$$$ AL_INI(3)=tath/sqrt(1+tath**2) c$$$ AL_INI(5)=(1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0. c$$$ c$$$c print*,'*******',AL_INI(5) c$$$ if(AL_INI(5).gt.defmax)goto 888 !next cloud c$$$ c$$$c print*,'alfaxz2, alfayz2 ' c$$$c $ ,alfaxz2_av(ixz),alfayz2_av(iyz) c$$$ c$$$* -------> INITIAL GUESS <------- c$$$c print*,'AL_INI ',(al_ini(i),i=1,5) c$$$ c$$$ if(DEBUG)then c$$$ print*,'1 >>> ',(cp_match(1,i),i=1,ncp_match(1)) c$$$ print*,'2 >>> ',(cp_match(2,i),i=1,ncp_match(2)) c$$$ print*,'3 >>> ',(cp_match(3,i),i=1,ncp_match(3)) c$$$ print*,'4 >>> ',(cp_match(4,i),i=1,ncp_match(4)) c$$$ print*,'5 >>> ',(cp_match(5,i),i=1,ncp_match(5)) c$$$ print*,'6 >>> ',(cp_match(6,i),i=1,ncp_match(6)) c$$$ endif c$$$ c$$$ do icp1=1,max(1,ncp_match(1)) c$$$ hit_plane(1)=icp1 c$$$ if(ncp_match(1).eq.0)hit_plane(1)=0 !-icp1 c$$$ c$$$ do icp2=1,max(1,ncp_match(2)) c$$$ hit_plane(2)=icp2 c$$$ if(ncp_match(2).eq.0)hit_plane(2)=0 !-icp2 c$$$ c$$$ do icp3=1,max(1,ncp_match(3)) c$$$ hit_plane(3)=icp3 c$$$ if(ncp_match(3).eq.0)hit_plane(3)=0 !-icp3 c$$$ c$$$ do icp4=1,max(1,ncp_match(4)) c$$$ hit_plane(4)=icp4 c$$$ if(ncp_match(4).eq.0)hit_plane(4)=0 !-icp4 c$$$ c$$$ do icp5=1,max(1,ncp_match(5)) c$$$ hit_plane(5)=icp5 c$$$ if(ncp_match(5).eq.0)hit_plane(5)=0 !-icp5 c$$$ c$$$ do icp6=1,max(1,ncp_match(6)) c$$$ hit_plane(6)=icp6 c$$$ if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6 c$$$ c$$$ c$$$ call track_init !init TRACK common c$$$ c$$$ do ip=1,nplanes !loop on planes c$$$ if(hit_plane(ip).ne.0)then c$$$ id=cp_match(ip,hit_plane(ip)) c$$$ is=is_cp(id) c$$$ icp=icp_cp(id) c$$$ if(ip_cp(id).ne.ip) c$$$ $ print*,'OKKIO!!' c$$$ $ ,'id ',id,is,icp c$$$ $ ,ip_cp(id),ip c$$$ icx=clx(ip,icp) c$$$ icy=cly(ip,icp) c$$$* ************************* c$$$ call xyz_PAM(icx,icy,is, c$$$ $ 'COG2','COG2',0.,0.) c$$$* ************************* c$$$* ----------------------------- c$$$ xgood(nplanes-ip+1)=1. c$$$ ygood(nplanes-ip+1)=1. c$$$ xm(nplanes-ip+1)=xPAM c$$$ ym(nplanes-ip+1)=yPAM c$$$ zm(nplanes-ip+1)=zPAM c$$$ resx(nplanes-ip+1)=resxPAM c$$$ resy(nplanes-ip+1)=resyPAM c$$$* ----------------------------- c$$$ endif c$$$ enddo !end loop on planes c$$$* ********************************************************** c$$$* ************************** FIT *** FIT *** FIT *** FIT *** c$$$* ********************************************************** c$$$ do i=1,5 c$$$ AL(i)=AL_INI(i) c$$$ enddo c$$$ ifail=0 !error flag in chi^2 computation c$$$ jstep=0 !number of minimization steps c$$$ call mini_2(jstep,ifail) c$$$ if(ifail.ne.0) then c$$$ if(DEBUG)then c$$$ print *, c$$$ $ '*** MINIMIZATION FAILURE *** ' c$$$ $ //'(mini_2 in clouds_to_ctrack)' c$$$ endif c$$$ chi2=-chi2 c$$$ endif c$$$* ********************************************************** c$$$* ************************** FIT *** FIT *** FIT *** FIT *** c$$$* ********************************************************** c$$$ c$$$ if(chi2.le.0.)goto 666 c$$$ c$$$* -------------------------- c$$$* STORE candidate TRACK INFO c$$$* -------------------------- c$$$ if(ntracks.eq.NTRACKSMAX)then c$$$ c$$$ if(DEBUG)print*, c$$$ $ '** warning ** number of candidate tracks '// c$$$ $ ' exceeds vector dimention ' c$$$ $ ,'( ',NTRACKSMAX,' )' c$$$c good2=.false. c$$$c goto 880 !fill ntp and go to next event c$$$ iflag=1 c$$$ return c$$$ endif c$$$ c$$$ ntracks = ntracks + 1 c$$$ c$$$c$$$ ndof=0 c$$$ do ip=1,nplanes c$$$c$$$ ndof=ndof c$$$c$$$ $ +int(xgood(ip)) c$$$c$$$ $ +int(ygood(ip)) c$$$ XV_STORE(ip,ntracks)=sngl(xv(ip)) c$$$ YV_STORE(ip,ntracks)=sngl(yv(ip)) c$$$ ZV_STORE(ip,ntracks)=sngl(zv(ip)) c$$$ XM_STORE(ip,ntracks)=sngl(xm(ip)) c$$$ YM_STORE(ip,ntracks)=sngl(ym(ip)) c$$$ ZM_STORE(ip,ntracks)=sngl(zm(ip)) c$$$ RESX_STORE(ip,ntracks)=sngl(resx(ip)) c$$$ RESY_STORE(ip,ntracks)=sngl(resy(ip)) c$$$ XV_STORE(ip,ntracks)=sngl(xv(ip)) c$$$ YV_STORE(ip,ntracks)=sngl(yv(ip)) c$$$ ZV_STORE(ip,ntracks)=sngl(zv(ip)) c$$$ AXV_STORE(ip,ntracks)=sngl(axv(ip)) c$$$ AYV_STORE(ip,ntracks)=sngl(ayv(ip)) c$$$ XGOOD_STORE(ip,ntracks)=sngl(xgood(ip)) c$$$ YGOOD_STORE(ip,ntracks)=sngl(ygood(ip)) c$$$ if(hit_plane(ip).ne.0)then c$$$ CP_STORE(nplanes-ip+1,ntracks)= c$$$ $ cp_match(ip,hit_plane(ip)) c$$$ else c$$$ CP_STORE(nplanes-ip+1,ntracks)=0 c$$$ endif c$$$ CLS_STORE(nplanes-ip+1,ntracks)=0 c$$$ do i=1,5 c$$$ AL_STORE(i,ntracks)=sngl(AL(i)) c$$$ enddo c$$$ enddo c$$$ c$$$c$$$ * Number of Degree Of Freedom c$$$c$$$ ndof=ndof-5 c$$$c$$$ * reduced chi^2 c$$$c$$$ rchi2=chi2/dble(ndof) c$$$ RCHI2_STORE(ntracks)=chi2 c$$$ c$$$* -------------------------------- c$$$* STORE candidate TRACK INFO - end c$$$* -------------------------------- c$$$ c$$$ 666 continue c$$$ enddo !end loop on cp in plane 6 c$$$ enddo !end loop on cp in plane 5 c$$$ enddo !end loop on cp in plane 4 c$$$ enddo !end loop on cp in plane 3 c$$$ enddo !end loop on cp in plane 2 c$$$ enddo !end loop on cp in plane 1 c$$$ c$$$ 888 continue c$$$ enddo !end loop on XZ couds c$$$ enddo !end loop on YZ couds c$$$ c$$$ if(ntracks.eq.0)then c$$$ iflag=1 c$$$ return c$$$ endif c$$$ c$$$ if(DEBUG)then c$$$ print*,'****** TRACK CANDIDATES ***********' c$$$ print*,'# R. chi2 RIG' c$$$ do i=1,ntracks c$$$ print*,i,' --- ',rchi2_store(i),' --- ' c$$$ $ ,1./abs(AL_STORE(5,i)) c$$$ enddo c$$$ print*,'***********************************' c$$$ endif c$$$ c$$$ c$$$ return c$$$ end c$$$ c$$$ c$$$*************************************************** c$$$* * c$$$* * c$$$* * c$$$* * c$$$* * c$$$* * c$$$************************************************** c$$$ c$$$ subroutine refine_track(ibest) c$$$ c$$$ include '../common/commontracker.f' c$$$ include '../common/common_momanhough.f' c$$$ include '../common/common_xyzPAM.f' c$$$ include '../common/common_mini_2.f' c$$$ include '../common/common_mech.f' c$$$ include '../common/momanhough_init.f' c$$$ include '../common/level1.f' c$$$ c$$$ logical DEBUG c$$$ common/dbg/DEBUG c$$$ c$$$* flag to chose PFA c$$$ character*10 PFA c$$$ common/FINALPFA/PFA c$$$ c$$$* ================================================= c$$$* new estimate of positions using ETA algorithm c$$$* and c$$$* search for new couples and single clusters to add c$$$* ================================================= c$$$ call track_init c$$$ do ip=1,nplanes !loop on planes c$$$ c$$$* ------------------------------------------------- c$$$* If the plane has been already included, it just c$$$* computes again the coordinates of the x-y couple c$$$* using improved PFAs c$$$* ------------------------------------------------- c$$$ if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and. c$$$ $ YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then c$$$ c$$$ id=CP_STORE(nplanes-ip+1,ibest) c$$$ c$$$ is=is_cp(id) c$$$ icp=icp_cp(id) c$$$ if(ip_cp(id).ne.ip) c$$$ $ print*,'OKKIO!!' c$$$ $ ,'id ',id,is,icp c$$$ $ ,ip_cp(id),ip c$$$ icx=clx(ip,icp) c$$$ icy=cly(ip,icp) c$$$ call xyz_PAM(icx,icy,is, c$$$c $ 'ETA2','ETA2', c$$$ $ PFA,PFA, c$$$ $ AXV_STORE(nplanes-ip+1,ibest), c$$$ $ AYV_STORE(nplanes-ip+1,ibest)) c$$$c$$$ call xyz_PAM(icx,icy,is, c$$$c$$$ $ 'COG2','COG2', c$$$c$$$ $ 0., c$$$c$$$ $ 0.) c$$$ xm(nplanes-ip+1) = xPAM c$$$ ym(nplanes-ip+1) = yPAM c$$$ zm(nplanes-ip+1) = zPAM c$$$ xgood(nplanes-ip+1) = 1 c$$$ ygood(nplanes-ip+1) = 1 c$$$ resx(nplanes-ip+1) = resxPAM c$$$ resy(nplanes-ip+1) = resyPAM c$$$ c$$$ dedxtrk(nplanes-ip+1) = (dedx(icx)+dedx(icy))/2. c$$$ c$$$* ------------------------------------------------- c$$$* If the plane has NOT been already included, c$$$* it tries to include a COUPLE or a single cluster c$$$* ------------------------------------------------- c$$$ else c$$$ c$$$ xgood(nplanes-ip+1)=0 c$$$ ygood(nplanes-ip+1)=0 c$$$ c$$$* -------------------------------------------------------------- c$$$* determine which ladder and sensor are intersected by the track c$$$ xP=XV_STORE(nplanes-ip+1,ibest) c$$$ yP=YV_STORE(nplanes-ip+1,ibest) c$$$ zP=ZV_STORE(nplanes-ip+1,ibest) c$$$ call whichsensor(ip,xP,yP,nldt,ist) c$$$* if the track hit the plane in a dead area, go to the next plane c$$$ if(nldt.eq.0.or.ist.eq.0)goto 133 c$$$* -------------------------------------------------------------- c$$$ c$$$ if(DEBUG)then c$$$ print*, c$$$ $ '------ Plane ',ip,' intersected on LADDER ',nldt c$$$ $ ,' SENSOR ',ist c$$$ print*, c$$$ $ '------ coord: ',XP,YP c$$$ endif c$$$ c$$$* =========================================== c$$$* STEP 1 >>>>>>> try to include a new couple c$$$* =========================================== c$$$c if(DEBUG)print*,'>>>> try to include a new couple' c$$$ distmin=1000000. c$$$ xmm = 0. c$$$ ymm = 0. c$$$ zmm = 0. c$$$ rxmm = 0. c$$$ rymm = 0. c$$$ idm = 0 !ID of the closer couple c$$$ distance=0. c$$$ do icp=1,ncp_plane(ip) !loop on couples on plane icp c$$$ icx=clx(ip,icp) c$$$ icy=cly(ip,icp) c$$$ if(LADDER(icx).ne.nldt.or. !If the ladder number does not match c$$$ $ cl_used(icx).eq.1.or. !or the X cluster is already used c$$$ $ cl_used(icy).eq.1.or. !or the Y cluster is already used c$$$ $ .false.)goto 1188 !then jump to next couple. c$$$* !Else compute the coordinates c$$$ call xyz_PAM(icx,icy,ist, c$$$ $ PFA,PFA, c$$$c $ 'ETA2','ETA2', c$$$ $ AXV_STORE(nplanes-ip+1,ibest), c$$$ $ AYV_STORE(nplanes-ip+1,ibest)) c$$$ c$$$ distance = distance_to(XP,YP) c$$$ id=id_cp(ip,icp,ist) c$$$ if(DEBUG)print*,'( couple ',id c$$$ $ ,' ) normalized distance ',distance c$$$ if(distance.lt.distmin)then c$$$ xmm = xPAM c$$$ ymm = yPAM c$$$ zmm = zPAM c$$$ rxmm = resxPAM c$$$ rymm = resyPAM c$$$ distmin = distance c$$$ idm = id c$$$ dedxmm = (dedx(icx)+dedx(icy))/2. c$$$ endif c$$$ 1188 continue c$$$ enddo !end loop on couples on plane icp c$$$ if(distmin.le.clinc)then c$$$* ----------------------------------- c$$$ xm(nplanes-ip+1) = xmm !<<< c$$$ ym(nplanes-ip+1) = ymm !<<< c$$$ zm(nplanes-ip+1) = zmm !<<< c$$$ xgood(nplanes-ip+1) = 1 !<<< c$$$ ygood(nplanes-ip+1) = 1 !<<< c$$$ resx(nplanes-ip+1)=rxmm !<<< c$$$ resy(nplanes-ip+1)=rymm !<<< c$$$ dedxtrk(nplanes-ip+1) = dedxmm !<<< c$$$* ----------------------------------- c$$$ CP_STORE(nplanes-ip+1,ibest)=idm c$$$ if(DEBUG)print*,'%%%% included couple ',idm c$$$ $ ,' (norm.dist.= ',distmin,', cut ',clinc,' )' c$$$ goto 133 !next plane c$$$ endif c$$$* ================================================ c$$$* STEP 2 >>>>>>> try to include a single cluster c$$$* either from a couple or single c$$$* ================================================ c$$$c if(DEBUG)print*,'>>>> try to include a new cluster' c$$$ distmin=1000000. c$$$ xmm_A = 0. !--------------------------- c$$$ ymm_A = 0. ! init variables that c$$$ zmm_A = 0. ! define the SINGLET c$$$ xmm_B = 0. ! c$$$ ymm_B = 0. ! c$$$ zmm_B = 0. ! c$$$ rxmm = 0. ! c$$$ rymm = 0. ! c$$$ iclm=0 !--------------------------- c$$$ distance=0. c$$$ c$$$*----- clusters inside couples ------------------------------------- c$$$ do icp=1,ncp_plane(ip) !loop on cluster inside couples c$$$ icx=clx(ip,icp) c$$$ icy=cly(ip,icp) c$$$ id=id_cp(ip,icp,ist) c$$$ if(LADDER(icx).ne.nldt)goto 11882 !if the ladder number does not match c$$$* !jump to the next couple c$$$*----- try cluster x ----------------------------------------------- c$$$ if(cl_used(icx).eq.1)goto 11881 !if the X cluster is already used c$$$* !jump to the Y cluster c$$$ call xyz_PAM(icx,0,ist, c$$$c $ 'ETA2','ETA2', c$$$ $ PFA,PFA, c$$$ $ AXV_STORE(nplanes-ip+1,ibest),0.) c$$$ distance = distance_to(XP,YP) c$$$c if(DEBUG)print*,'normalized distance ',distance c$$$ if(DEBUG)print*,'( cl-X ',icx c$$$ $ ,' in cp ',id,' ) normalized distance ',distance c$$$ if(distance.lt.distmin)then c$$$ xmm_A = xPAM_A c$$$ ymm_A = yPAM_A c$$$ zmm_A = zPAM_A c$$$ xmm_B = xPAM_B c$$$ ymm_B = yPAM_B c$$$ zmm_B = zPAM_B c$$$ rxmm = resxPAM c$$$ rymm = resyPAM c$$$ distmin = distance c$$$ iclm = icx c$$$ dedxmm = dedx(icx) c$$$ endif c$$$11881 continue c$$$*----- try cluster y ----------------------------------------------- c$$$ if(cl_used(icy).eq.1)goto 11882 !if the Y cluster is already used c$$$* !jump to the next couple c$$$ call xyz_PAM(0,icy,ist, c$$$c $ 'ETA2','ETA2', c$$$ $ PFA,PFA, c$$$ $ 0.,AYV_STORE(nplanes-ip+1,ibest)) c$$$ distance = distance_to(XP,YP) c$$$ if(DEBUG)print*,'( cl-Y ',icy c$$$ $ ,' in cp ',id,' ) normalized distance ',distance c$$$ if(distance.lt.distmin)then c$$$ xmm_A = xPAM_A c$$$ ymm_A = yPAM_A c$$$ zmm_A = zPAM_A c$$$ xmm_B = xPAM_B c$$$ ymm_B = yPAM_B c$$$ zmm_B = zPAM_B c$$$ rxmm = resxPAM c$$$ rymm = resyPAM c$$$ distmin = distance c$$$ iclm = icy c$$$ dedxmm = dedx(icy) c$$$ endif c$$$11882 continue c$$$ enddo !end loop on cluster inside couples c$$$*----- single clusters ----------------------------------------------- c$$$ do ic=1,ncls(ip) !loop on single clusters c$$$ icl=cls(ip,ic) c$$$ if(cl_used(icl).eq.1.or. !if the cluster is already used c$$$ $ LADDER(icl).ne.nldt.or. !or the ladder number does not match c$$$ $ .false.)goto 18882 !jump to the next singlet c$$$ if(mod(VIEW(icl),2).eq.0)then!<---- X view c$$$ call xyz_PAM(icl,0,ist, c$$$c $ 'ETA2','ETA2', c$$$ $ PFA,PFA, c$$$ $ AXV_STORE(nplanes-ip+1,ibest),0.) c$$$ else !<---- Y view c$$$ call xyz_PAM(0,icl,ist, c$$$c $ 'ETA2','ETA2', c$$$ $ PFA,PFA, c$$$ $ 0.,AYV_STORE(nplanes-ip+1,ibest)) c$$$ endif c$$$ c$$$ distance = distance_to(XP,YP) c$$$ if(DEBUG)print*,'( cl-s ',icl c$$$ $ ,' ) normalized distance ',distance c$$$ if(distance.lt.distmin)then c$$$ xmm_A = xPAM_A c$$$ ymm_A = yPAM_A c$$$ zmm_A = zPAM_A c$$$ xmm_B = xPAM_B c$$$ ymm_B = yPAM_B c$$$ zmm_B = zPAM_B c$$$ rxmm = resxPAM c$$$ rymm = resyPAM c$$$ distmin = distance c$$$ iclm = icl c$$$ dedxmm = dedx(icl) c$$$ endif c$$$18882 continue c$$$ enddo !end loop on single clusters c$$$ c$$$ if(distmin.le.clinc)then c$$$ c$$$ CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<< c$$$* ---------------------------- c$$$ if(mod(VIEW(iclm),2).eq.0)then c$$$ XGOOD(nplanes-ip+1)=1. c$$$ resx(nplanes-ip+1)=rxmm c$$$ if(DEBUG)print*,'%%%% included X-cl ',iclm c$$$ $ ,' ( norm.dist.= ',distmin,', cut ',clinc,' )' c$$$ else c$$$ YGOOD(nplanes-ip+1)=1. c$$$ resy(nplanes-ip+1)=rymm c$$$ if(DEBUG)print*,'%%%% included Y-cl ',iclm c$$$ $ ,' ( norm.dist.= ',distmin,', cut ',clinc,' )' c$$$ endif c$$$* ---------------------------- c$$$ xm_A(nplanes-ip+1) = xmm_A c$$$ ym_A(nplanes-ip+1) = ymm_A c$$$ xm_B(nplanes-ip+1) = xmm_B c$$$ ym_B(nplanes-ip+1) = ymm_B c$$$ zm(nplanes-ip+1) = (zmm_A+zmm_B)/2. c$$$ dedxtrk(nplanes-ip+1) = dedxmm !<<< c$$$* ---------------------------- c$$$ endif c$$$ endif c$$$ 133 continue c$$$ enddo !end loop on planes c$$$ c$$$ c$$$ c$$$ return c$$$ end c$$$ c$$$*************************************************** c$$$* * c$$$* * c$$$* * c$$$* * c$$$* * c$$$* * c$$$************************************************** c$$$ c$$$ subroutine clean_XYclouds(ibest,iflag) c$$$ c$$$ include '../common/commontracker.f' c$$$ include '../common/common_momanhough.f' c$$$ include '../common/momanhough_init.f' c$$$c include '../common/calib.f' c$$$c include '../common/level1.f' c$$$ c$$$ logical DEBUG c$$$ common/dbg/DEBUG c$$$ c$$$ c$$$ do ip=1,nplanes !loop on planes c$$$ c$$$ id=CP_STORE(nplanes-ip+1,ibest) c$$$ icl=CLS_STORE(nplanes-ip+1,ibest) c$$$ if(id.ne.0.or.icl.ne.0)then c$$$ if(id.ne.0)then c$$$ iclx=clx(ip,icp_cp(id)) c$$$ icly=cly(ip,icp_cp(id)) c$$$ cl_used(iclx)=1 !tag used clusters c$$$ cl_used(icly)=1 !tag used clusters c$$$ elseif(icl.ne.0)then c$$$ cl_used(icl)=1 !tag used clusters c$$$ endif c$$$ c$$$c if(DEBUG)then c$$$c print*,ip,' <<< ',id c$$$c endif c$$$* ----------------------------- c$$$* remove the couple from clouds c$$$* remove also vitual couples containing the c$$$* selected clusters c$$$* ----------------------------- c$$$ do icp=1,ncp_plane(ip) c$$$ if( c$$$ $ clx(ip,icp).eq.iclx c$$$ $ .or. c$$$ $ clx(ip,icp).eq.icl c$$$ $ .or. c$$$ $ cly(ip,icp).eq.icly c$$$ $ .or. c$$$ $ cly(ip,icp).eq.icl c$$$ $ )then c$$$ id=id_cp(ip,icp,1) c$$$ if(DEBUG)then c$$$ print*,ip,' <<< cp ',id c$$$ $ ,' ( cl-x ' c$$$ $ ,clx(ip,icp) c$$$ $ ,' cl-y ' c$$$ $ ,cly(ip,icp),' ) --> removed' c$$$ endif c$$$* ----------------------------- c$$$* remove the couple from clouds c$$$ do iyz=1,nclouds_yz c$$$ if(cpcloud_yz(iyz,abs(id)).ne.0)then c$$$ ptcloud_yz(iyz)=ptcloud_yz(iyz)-1 c$$$ cpcloud_yz(iyz,abs(id))=0 c$$$ endif c$$$ enddo c$$$ do ixz=1,nclouds_xz c$$$ if(cpcloud_xz(ixz,abs(id)).ne.0)then c$$$ ptcloud_xz(ixz)=ptcloud_xz(ixz)-1 c$$$ cpcloud_xz(ixz,abs(id))=0 c$$$ endif c$$$ enddo c$$$* ----------------------------- c$$$ endif c$$$ enddo c$$$ c$$$ endif c$$$ enddo !end loop on planes c$$$ c$$$ return c$$$ end c$$$ c$$$ c$$$ c$$$ c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** c$$$ real function fbad_cog(ncog,ic) c$$$ c$$$ c$$$ include '../common/commontracker.f' c$$$ include '../common/level1.f' c$$$ include '../common/calib.f' c$$$ c$$$* --> signal of the central strip c$$$ sc = CLSIGNAL(INDMAX(ic)) !center c$$$ c$$$* signal of adjacent strips c$$$* --> left c$$$ sl1 = 0 !left 1 c$$$ if( c$$$ $ (INDMAX(ic)-1).ge.INDSTART(ic) c$$$ $ ) c$$$ $ sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1)) c$$$ c$$$ sl2 = 0 !left 2 c$$$ if( c$$$ $ (INDMAX(ic)-2).ge.INDSTART(ic) c$$$ $ ) c$$$ $ sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2)) c$$$ c$$$* --> right c$$$ sr1 = 0 !right 1 c$$$ if( c$$$ $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1)) c$$$ $ .or. c$$$ $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH) c$$$ $ ) c$$$ $ sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1)) c$$$ c$$$ sr2 = 0 !right 2 c$$$ if( c$$$ $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1)) c$$$ $ .or. c$$$ $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH) c$$$ $ ) c$$$ $ sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2)) c$$$ c$$$ c$$$ if(mod(int(VIEW(ic)),2).eq.1)then !Y-view c$$$ f = 4. c$$$ si = 8.4 c$$$ else !X-view c$$$ f = 6. c$$$ si = 3.9 c$$$ endif c$$$ c$$$ fbad_cog = 1. c$$$ f0 = 1 c$$$ f1 = 1 c$$$ f2 = 1 c$$$ f3 = 1 c$$$ if(sl1.gt.sr1.and.sl1.gt.0.)then c$$$ c$$$ if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)) ).eq.0)f0=f c$$$ if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f c$$$c if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f c$$$ c$$$ if(ncog.eq.2.and.sl1.ne.0)then c$$$ fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.) c$$$ elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then c$$$ fbad_cog = 1. c$$$ elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then c$$$ fbad_cog = 1. c$$$ else c$$$ fbad_cog = 1. c$$$ endif c$$$ c$$$ elseif(sl1.le.sr1.and.sr1.gt.0.)then c$$$ c$$$ c$$$ if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)) ).eq.0)f0=f c$$$ if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f c$$$c if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f c$$$ c$$$ if(ncog.eq.2.and.sr1.ne.0)then c$$$ fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.) c$$$ elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then c$$$ fbad_cog = 1. c$$$ elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then c$$$ fbad_cog = 1. c$$$ else c$$$ fbad_cog = 1. c$$$ endif c$$$ c$$$ endif c$$$ c$$$ fbad_cog = sqrt(fbad_cog) c$$$ c$$$ return c$$$ end