--- DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2006/11/08 16:42:28 1.13 +++ DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2006/11/21 14:00:40 1.14 @@ -18,7 +18,6 @@ include 'common_xyzPAM.f' include 'common_mini_2.f' include 'calib.f' -c include 'level1.f' include 'level2.f' c include 'momanhough_init.f' @@ -423,178 +422,6 @@ end - - -c$$$************************************************************ -c$$$ -c$$$ subroutine readmipparam -c$$$ -c$$$ include 'commontracker.f' -c$$$ include 'calib.f' -c$$$ -c$$$ character*60 fname_param -c$$$ 201 format('trk-LADDER',i1,'-mip.dat') -c$$$ do ilad=1,nladders_view -c$$$ write(fname_param,201)ilad -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*,'READMIPPARAM: *** Error in opening file ***' -c$$$ return -c$$$ endif -c$$$ do iv=1,nviews -c$$$ read(10,* -c$$$ $ ,IOSTAT=iostat -c$$$ $ )pip, -c$$$ $ mip(int(pip),ilad) -c$$$c print*,ilad,iv,pip,mip(int(pip),ilad) -c$$$ enddo -c$$$ close(10) -c$$$ enddo -c$$$ -c$$$ return -c$$$ end -c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** -c$$$ subroutine readchargeparam -c$$$ -c$$$ -c$$$ include 'commontracker.f' -c$$$ include 'calib.f' -c$$$ -c$$$ character*60 fname_param -c$$$ 201 format('charge-l',i1,'.dat') -c$$$ do ilad=1,nladders_view -c$$$ write(fname_param,201)ilad -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$$$ enddo -c$$$ -c$$$ return -c$$$ end -c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** -c$$$ subroutine readetaparam -c$$$* ----------------------------------------- -c$$$* read eta2,3,4 calibration parameters -c$$$* and fill variables: -c$$$* -c$$$* eta2(netabin,nladders_view,nviews) -c$$$* eta3(2*netabin,nladders_view,nviews) -c$$$* eta4(2*netabin,nladders_view,nviews) -c$$$* -c$$$ include 'commontracker.f' -c$$$ include 'calib.f' -c$$$ -c$$$ character*40 fname_binning -c$$$ character*40 fname_param -c$$$c character*120 cmd1 -c$$$c character*120 cmd2 -c$$$ -c$$$ -c$$$******retrieve ANGULAR BINNING info -c$$$ fname_binning='binning.dat' -c$$$ print *,'Opening file: ',fname_binning -c$$$ open(10, -c$$$ $ FILE='./bin-aux/'//fname_binning(1:LNBLNK(fname_binning)) -c$$$ $ ,STATUS='UNKNOWN' -c$$$ $ ,IOSTAT=iostat -c$$$ $ ) -c$$$ if(iostat.ne.0)then -c$$$ print*,'READETAPARAM: *** Error in opening file ***' -c$$$ return -c$$$ endif -c$$$ print*,'---- ANGULAR BINNING ----' -c$$$ print*,'Bin - angL - angR' -c$$$ 101 format(i2,' ',f6.2,' ',f6.2) -c$$$ do ibin=1,nangmax -c$$$ read(10,* -c$$$ $ ,IOSTAT=iostat -c$$$ $ )xnn,angL(ibin),angR(ibin) -c$$$ if(iostat.ne.0)goto 1000 -c$$$ write(*,101)int(xnn),angL(ibin),angR(ibin) -c$$$ enddo -c$$$ 1000 nangbin=int(xnn) -c$$$ close(10) -c$$$ print*,'-------------------------' -c$$$ -c$$$ -c$$$ -c$$$ do ieta=2,4 !loop on eta 2,3,4 -c$$$******retrieve correction parameters -c$$$ 200 format(' Opening eta',i1,' files...') -c$$$ write(*,200)ieta -c$$$ -c$$$ 201 format('eta',i1,'-bin',i1,'-l',i1,'.dat') -c$$$ 202 format('eta',i1,'-bin',i2,'-l',i1,'.dat') -c$$$ do iang=1,nangbin -c$$$ do ilad=1,nladders_view -c$$$ if(iang.lt.10)write(fname_param,201)ieta,iang,ilad -c$$$ if(iang.ge.10)write(fname_param,202)ieta,iang,ilad -c$$$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*,'READETAPARAM: *** Error in opening file ***' -c$$$ return -c$$$ endif -c$$$ do ival=1,netavalmax -c$$$ if(ieta.eq.2)read(10,* -c$$$ $ ,IOSTAT=iostat -c$$$ $ ) -c$$$ $ eta2(ival,iang), -c$$$ $ (feta2(ival,iv,ilad,iang),iv=1,nviews) -c$$$ if(ieta.eq.3)read(10,* -c$$$ $ ,IOSTAT=iostat -c$$$ $ ) -c$$$ $ eta3(ival,iang), -c$$$ $ (feta3(ival,iv,ilad,iang),iv=1,nviews) -c$$$ if(ieta.eq.4)read(10,* -c$$$ $ ,IOSTAT=iostat -c$$$ $ ) -c$$$ $ eta4(ival,iang), -c$$$ $ (feta4(ival,iv,ilad,iang),iv=1,nviews) -c$$$ if(iostat.ne.0)then -c$$$ netaval=ival-1 -c$$$c$$$ if(eta2(1,iang).ne.-eta2(netaval,iang)) -c$$$c$$$ $ print*,'**** ERROR on parameters !!! ****' -c$$$ goto 2000 -c$$$ endif -c$$$ enddo -c$$$ 2000 close(10) -c$$$* print*,'... done' -c$$$ enddo -c$$$ enddo -c$$$ -c$$$ enddo !end loop on eta 2,3,4 -c$$$ -c$$$ -c$$$ return -c$$$ end -c$$$ - ************************************************************ ************************************************************ @@ -1444,274 +1271,6 @@ ************************************************************************* ************************************************************************* ************************************************************************* -c$$$ subroutine book_debug -c$$$ -c$$$ include 'commontracker.f' -c$$$ include 'common_momanhough.f' -c$$$ include 'common_level2debug.f' -c$$$ -c$$$ character*35 block1,block2,block3!,block4 -c$$$ $ ,block5!,block6 -c$$$ -c$$$* * * * * * * * * * * * * * * * * * * * * * * * -c$$$* HOUGH TRANSFORM PARAMETERS -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$$$ $ ,DB_CLOUD(NDBLT):I -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$$$ $ ,TR_CLOUD(NTRPT):I -c$$$ $ ') -c$$$ -c$$$ -c$$$ 413 format('NCLOUDS_YZ:I::[0,',I4,']') -c$$$c$$$ 414 format('DB_CLOUD(',I4,'):I') -c$$$ write(block3,413) ncloyz_max -c$$$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$$$c$$$ $ ,'//block4 -c$$$ $ ) -c$$$ -c$$$ 415 format('NCLOUDS_XZ:I::[0,',I4,']') -c$$$c$$$ 416 format('TR_CLOUD(',I5,'):I') -c$$$ write(block5,415) ncloxz_max -c$$$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$$$c$$$ $ ,'//block6 -c$$$ $ ) -c$$$ -c$$$ -c$$$ return -c$$$ end -***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...*** -* -* -* -* -* -* -* -* -* -***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...*** -c$$$ subroutine book_level2 -c$$$c***************************************************** -c$$$cccccc 11/9/2005 modified by david fedele -c$$$cccccc 07/10/2005 modified by elena vannuccini --> (2) -c$$$c***************************************************** -c$$$ -c$$$ include 'commontracker.f' -c$$$ include 'common_momanhough.f' -c$$$ include '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$$$c***************************************************** -c$$$cccccc 11/9/2005 modified by david fedele -c$$$c call HBNAME(ntp_level2,'EVENT',good2,'GOOD2:L,NEV2:I') -c$$$cccccc 06/10/2005 modified by elena vannuccini -c$$$c call HBNAME(ntp_level2,'GENERAL',good2,'GOOD2:L,NEV2:I -c$$$c $ ,WHIC_CALIB:I,SWCODE:I') -c$$$ call HBNAME(ntp_level2,'GENERAL',good2,'GOOD2:L,NEV2:I -c$$$ $ ,WHICH_CALIB:I,SWCODE:I,CRC(12):L') -c$$$c********************************************************* -c$$$ -c$$$ -c$$$c$$$# ifndef TEST2003 -c$$$c$$$ -c$$$c$$$ call HBNAME(ntp_level2,'CPU',pkt_type -c$$$c$$$ $ ,'PKT_TYPE:I::[0,50] -c$$$c$$$ $ ,PKT_NUM:I -c$$$c$$$ $ ,OBT:I'// -c$$$c$$$c******************************************************** -c$$$c$$$cccccc 11/9/2005 modified by david fedele -c$$$c$$$c $ ,WHICH_CALIB:I::[0,50]') -c$$$c$$$ $ ',CPU_CRC:L') -c$$$c$$$c******************************************************** -c$$$c$$$ -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$$$c***************************************************** -c$$$cccccc 11/9/2005 modified by david fedele -c$$$c $ ,DEDXP(6,NTRK):R'// -c$$$c $ ') -c$$$ $ ',DEDX_X(6,NTRK):R -c$$$ $ ,DEDX_Y(6,NTRK):R'// -c$$$c**************************************************** -c$$$cccccc 06/10/2005 modified by elena vannuccini -c$$$c $ ,CRC(12):L -c$$$ $ ',BdL(NTRK):R' -c$$$ $ ) -c$$$c**************************************************** -c$$$ -c$$$ -c$$$ call HBNAME(ntp_level2,'SINGLETX',nclsx, -c$$$c***************************************************** -c$$$cccccc 11/9/2005 modified by david fedele -c$$$c $ 'NCLSX(6):I,NCLSY(6):I') -c$$$ $ 'NCLSX:I::[0,500],PLANEX(NCLSX):I -c$$$ $ ,XS(2,NCLSX):R,SGNLXS(NCLSX):R') !(2) -c$$$c $ ,XS(NCLSX):R,SGNLXS(NCLSX):R') !(2) -c$$$ call HBNAME(ntp_level2,'SINGLETY',nclsy, -c$$$ $ 'NCLSY:I::[0,500],PLANEY(NCLSY):I -c$$$ $ ,YS(2,NCLSY):R,SGNLYS(NCLSY):R') !(2) -c$$$c $ ,YS(NCLSY):R,SGNLYS(NCLSY):R') !(2) -c$$$ return -c$$$ end - -c$$$ subroutine fill_level2_clouds -c$$$c***************************************************** -c$$$c 29/11/2005 created by elena vannuccini -c$$$c***************************************************** -c$$$ -c$$$* ------------------------------------------------------- -c$$$* This routine fills the variables related to the hough -c$$$* transform, for the debig n-tuple -c$$$* ------------------------------------------------------- -c$$$ -c$$$ include 'commontracker.f' -c$$$ include 'common_momanhough.f' -c$$$ include 'common_level2debug.f' -c$$$ include 'level2.f' -c$$$ -c$$$ good2_nt=.true.!good2 -c$$$c nev2_nt=nev2 -c$$$ -c$$$ if(.false. -c$$$ $ .or.ntrpt.gt.ntrpt_max_nt -c$$$ $ .or.ndblt.gt.ndblt_max_nt -c$$$ $ .or.NCLOUDS_XZ.gt.ncloxz_max -c$$$ $ .or.NCLOUDS_yZ.gt.ncloyz_max -c$$$ $ )then -c$$$ good2_nt=.false. -c$$$ ntrpt_nt=0 -c$$$ ndblt_nt=0 -c$$$ NCLOUDS_XZ_nt=0 -c$$$ NCLOUDS_YZ_nt=0 -c$$$ else -c$$$ ndblt_nt=ndblt -c$$$ ntrpt_nt=ntrpt -c$$$ if(ndblt.ne.0)then -c$$$ do id=1,ndblt -c$$$ alfayz1_nt(id)=alfayz1(id) !Y0 -c$$$ alfayz2_nt(id)=alfayz2(id) !tg theta-yz -c$$$c db_cloud_nt(id)=db_cloud(id) -c$$$ enddo -c$$$ endif -c$$$ if(ndblt.ne.0)then -c$$$ do it=1,ntrpt -c$$$ alfaxz1_nt(it)=alfaxz1(it) !X0 -c$$$ alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz -c$$$ alfaxz3_nt(it)=alfaxz3(it) !1/r -c$$$c tr_cloud_nt(it)=tr_cloud(it) -c$$$ enddo -c$$$ endif -c$$$ nclouds_yz_nt=nclouds_yz -c$$$ nclouds_xz_nt=nclouds_xz -c$$$ if(nclouds_yz.ne.0)then -c$$$ nnn=0 -c$$$ do iyz=1,nclouds_yz -c$$$ ptcloud_yz_nt(iyz)=ptcloud_yz(iyz) -c$$$ alfayz1_av_nt(iyz)=alfayz1_av(iyz) -c$$$ alfayz2_av_nt(iyz)=alfayz2_av(iyz) -c$$$ nnn=nnn+ptcloud_yz(iyz) -c$$$ enddo -c$$$ do ipt=1,nnn -c$$$ db_cloud_nt(ipt)=db_cloud(ipt) -c$$$ enddo -c$$$c print*,'#### ntupla #### ptcloud_yz ' -c$$$c $ ,(ptcloud_yz(i),i=1,nclouds_yz) -c$$$c print*,'#### ntupla #### db_cloud ' -c$$$c $ ,(db_cloud(i),i=1,nnn) -c$$$ endif -c$$$ if(nclouds_xz.ne.0)then -c$$$ nnn=0 -c$$$ do ixz=1,nclouds_xz -c$$$ ptcloud_xz_nt(ixz)=ptcloud_xz(ixz) -c$$$ alfaxz1_av_nt(ixz)=alfaxz1_av(ixz) -c$$$ alfaxz2_av_nt(ixz)=alfaxz2_av(ixz) -c$$$ alfaxz3_av_nt(ixz)=alfaxz3_av(ixz) -c$$$ nnn=nnn+ptcloud_xz(ixz) -c$$$ enddo -c$$$ do ipt=1,nnn -c$$$ tr_cloud_nt(ipt)=tr_cloud(ipt) -c$$$ enddo -c$$$c print*,'#### ntupla #### ptcloud_xz ' -c$$$c $ ,(ptcloud_xz(i),i=1,nclouds_xz) -c$$$c print*,'#### ntupla #### tr_cloud ' -c$$$c $ ,(tr_cloud(i),i=1,nnn) -c$$$ endif -c$$$ endif -c$$$ end *************************************************** @@ -1769,10 +1328,10 @@ enddo * mask views with too many clusters do iv=1,nviews - if( ncl_view(iv).gt. nclustermax)then + if( ncl_view(iv).gt. nclusterlimit)then mask_view(iv) = 1 if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > ' - $ ,nclustermax,' on view ', iv,' --> masked!' + $ ,nclusterlimit,' on view ', iv,' --> masked!' endif enddo @@ -1910,41 +1469,25 @@ endif * ------------------> COUPLE <------------------ -* 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 - + ncp_plane(nplx) = ncp_plane(nplx) + 1 + clx(nplx,ncp_plane(nplx))=icx + cly(nply,ncp_plane(nplx))=icy + cl_single(icx)=0 + cl_single(icy)=0 + if(ncp_plane(nplx).eq.ncouplemax)then if(verbose)print*, $ '** warning ** number of identified '// $ 'couples on plane ',nplx, $ 'exceeds vector dimention ' - $ ,'( ',ncouplemax,' )' -c good2=.false. -c goto 880 !fill ntp and go to next event + $ ,'( ',ncouplemax,' ) --> masked!' mask_view(nviewx(nplx)) = 2 mask_view(nviewy(nply)) = 2 -c iflag=1 -c return endif - - ncp_plane(nplx) = ncp_plane(nplx) + 1 - clx(nplx,ncp_plane(nplx))=icx - cly(nply,ncp_plane(nplx))=icy - cl_single(icx)=0 - cl_single(icy)=0 - endif * ---------------------------------------------- + endif + 20 continue enddo !end loop on clusters(Y) @@ -1971,16 +1514,6 @@ do ip=1,6 ncp_tot = ncp_tot + ncp_plane(ip) enddo -c if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!) - -c$$$ if(ncp_tot.gt.ncp_max)then -c$$$ if(verbose)print*, -c$$$ $ '** warning ** number of identified '// -c$$$ $ 'couples exceeds upper limit for Hough tr. ' -c$$$ $ ,'( ',ncp_max,' )' -c$$$ iflag=1 -c$$$ return -c$$$ endif return end @@ -1993,240 +1526,15 @@ * * * * ************************************************** -c$$$ subroutine cl_to_couples_nocharge(iflag) -c$$$ -c$$$ include 'commontracker.f' -c$$$ include 'level1.f' -c$$$ include 'common_momanhough.f' -c$$$c include 'momanhough_init.f' -c$$$ include 'calib.f' -c$$$c include 'level1.f' -c$$$ -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$$$ if(badcl.eq.0)then !<<<<<<<<<<<<<< BAD cut -c$$$ cl_single(icx)=0 !<<<<<<<<<<<<<< BAD cut -c$$$ goto 10 !<<<<<<<<<<<<<< BAD cut -c$$$ endif !<<<<<<<<<<<<<< BAD cut -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$$$ if(badcl.eq.0)then !<<<<<<<<<<<<<< BAD cut -c$$$ cl_single(icy)=0 !<<<<<<<<<<<<<< BAD cut -c$$$ goto 20 !<<<<<<<<<<<<<< BAD cut -c$$$ endif !<<<<<<<<<<<<<< BAD cut -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$$$* =========================================================== -c$$$* this version of the subroutine is used for the calibration -c$$$* thus charge-correlation selection is obviously removed -c$$$* =========================================================== -c$$$c$$$ ddd=(dedx(icy) -c$$$c$$$ $ -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx)) -c$$$c$$$ ddd=ddd/sqrt(kch(nplx,nldx)**2+1) -c$$$c$$$ cut=chcut*sch(nplx,nldx) -c$$$c$$$ if(abs(ddd).gt.cut)goto 20 !charge not consistent -c$$$* =========================================================== -c$$$ -c$$$ -c$$$* ------------------> COUPLE <------------------ -c$$$* check to do not overflow vector dimentions -c$$$c$$$ if(ncp_plane(nplx).gt.ncouplemax)then -c$$$c$$$ if(DEBUG)print*, -c$$$c$$$ $ ' ** warning ** number of identified'// -c$$$c$$$ $ ' couples on plane ',nplx, -c$$$c$$$ $ ' exceeds vector dimention'// -c$$$c$$$ $ ' ( ',ncouplemax,' )' -c$$$c$$$c good2=.false. -c$$$c$$$c goto 880 !fill ntp and go to next event -c$$$c$$$ iflag=1 -c$$$c$$$ return -c$$$c$$$ endif -c$$$ -c$$$ if(ncp_plane(nplx).eq.ncouplemax)then -c$$$ if(verbose)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(verbose)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$$$ - -*************************************************** -* * -* * -* * -* * -* * -* * -************************************************** subroutine cp_to_doubtrip(iflag) -c***************************************************** -c 02/02/2006 modified by Elena Vannuccini --> (1) -c***************************************************** include 'commontracker.f' include 'level1.f' include 'common_momanhough.f' -c include 'momanhough_init.f' include 'common_xyzPAM.f' include 'common_mini_2.f' include 'calib.f' -c include 'level1.f' * output flag @@ -2260,13 +1568,26 @@ * ----------------------------- +* -------------------------------------------- +* put a limit to the maximum number of couples +* per plane, in order to apply hough transform +* (couples recovered during track refinement) +* -------------------------------------------- + do ip=1,nplanes + if(ncp_plane(ip).gt.ncouplelimit)then + mask_view(nviewx(ip)) = 8 + mask_view(nviewy(ip)) = 8 + endif + enddo + ndblt=0 !number of doublets ntrpt=0 !number of triplets do ip1=1,(nplanes-1) !loop on planes - COPPIA 1 - do is1=1,2 !loop on sensors - COPPIA 1 - + if( mask_view(nviewx(ip1)).ne.0 .or. + $ mask_view(nviewy(ip1)).ne.0 )goto 10 !skip plane + do is1=1,2 !loop on sensors - COPPIA 1 do icp1=1,ncp_plane(ip1) !loop on COPPIA 1 icx1=clx(ip1,icp1) icy1=cly(ip1,icp1) @@ -2276,7 +1597,10 @@ ym1=yPAM zm1=zPAM c print*,'***',is1,xm1,ym1,zm1 + do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2 + if( mask_view(nviewx(ip2)).ne.0 .or. + $ mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane do is2=1,2 !loop on sensors -ndblt COPPIA 2 do icp2=1,ncp_plane(ip2) !loop on COPPIA 2 @@ -2337,8 +1661,11 @@ * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - if(ip2.eq.nplanes)goto 30 !no possible combination with 3 couples + if(ip2.eq.nplanes)goto 31 !no possible combination with 3 couples + do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3 + if( mask_view(nviewx(ip3)).ne.0 .or. + $ mask_view(nviewy(ip3)).ne.0 )goto 30 !skip plane do is3=1,2 !loop on sensors - COPPIA 3 do icp3=1,ncp_plane(ip3) !loop on COPPIA 3 @@ -2415,15 +1742,18 @@ endif enddo !end loop on COPPIA 3 enddo !end loop on sensors - COPPIA 3 + 30 continue enddo !end loop on planes - COPPIA 3 - 30 continue + 31 continue 1 enddo !end loop on COPPIA 2 enddo !end loop on sensors - COPPIA 2 + 20 continue enddo !end loop on planes - COPPIA 2 enddo !end loop on COPPIA1 enddo !end loop on sensors - COPPIA 1 + 10 continue enddo !end loop on planes - COPPIA 1 if(DEBUG)then @@ -3684,103 +3014,6 @@ -c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** -c$$$ real function fbad_cog(ncog,ic) -c$$$ -c$$$ -c$$$ include 'commontracker.f' -c$$$ include 'level1.f' -c$$$ include '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 -c$$$ - * **************************************************** @@ -3865,7 +3098,7 @@ alfayz1_nt(idb)=0 alfayz2_nt(idb)=0 enddo - do itr=1,ntrpl_max_nt + do itr=1,ntrpt_max_nt tr_cloud_nt(itr)=0 alfaxz1_nt(itr)=0 alfaxz2_nt(itr)=0 @@ -3894,7 +3127,7 @@ alfayz1(idb)=0 alfayz2(idb)=0 enddo - do itr=1,ntrpl_max + do itr=1,ntrpt_max tr_cloud(itr)=0 cpxz1(itr)=0 cpxz2(itr)=0