C =============================================== subroutine MUONFIND(nomu,mupt,muth,muph,NrMips) C =============================================== c ====================================================================== c 22.10.97 : Find calorimeter MIPS with matching tracks for ZES. c D.Waters (Oxford) c 23.11.97 : Modified to store nCalMip in ZES c A.Quadt (Oxford) c 30.11.98 : Modified to use VCTDCA rather than TRKCLU. c ====================================================================== c 1 2 3 4 5 6 7 implicit none c Calorimeter MIP related sequences : c ----------------------------------- #include "islandc.inc" #include "calmipc.inc" c For vertex position and tracks : c -------------------------------- #include "partap.inc" #include "vctvtx.inc" #include "vctpar.inc" #include "vctrhl.inc" c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c Output from subroutine MUONFIND : c ================================= c Number of muon candidates passing muon selection cuts : c ------------------------------------------------------- integer nomu c Pt of the 2 highest Pt candidates : c ----------------------------------- real mupt(2) c Polar angle of the 2 highest Pt candidates : c -------------------------------------------- real muth(2) c Azimuthal angle of the 2 highest Pt candidates : c ------------------------------------------------ real muph(2) c Nr of mips deposited in the calorimeter : c ----------------------------------------- integer NrMips c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c Muon selection cuts : c ===================== c Select muon with Pt above `pt_cut' GeV : c ---------------------------------------- real pt_cut parameter (pt_cut=2.0) c DCA (closest track) < `dca' cm. : c --------------------------------- real dca parameter (dca=30.) c Isolation parameter (veto if other tracks have DCA < `iso' cm. : c ---------------------------------------------------------------- real iso parameter (iso=0.) c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ integer i,j,mipno,FMode,Mode,Nfound,iErr,indx(200),num_match integer sort_mode,sort_nway,nsort,vctrhl_index(200) real dist_app(200),vertex(3) real trk_theta,pt_sort(10),th_sort(10),ph_sort(10) logical ok c For calls to VCTDCA : c --------------------- integer vcid real calpos(3),calerr(3),vctdca_dca,covprojtrk,covprojcal c ---------------------------------------------------------------------- c Initialise some variables : c --------------------------- nomu=0 do i=1,2 mupt(i)=0. muth(i)=0. muph(i)=0. enddo do i=1,10 pt_sort(i)=0. th_sort(i)=0. ph_sort(i)=0. enddo c GLOMU island, archipelago and MIP finding : c =========================================== call Zeroit call CellFill call Island call IslFill c Set error flag to zero first : iCalMipErr=0 call CalMip if (iCalMipErr.NE.0) then c write(6,*) 'Error MUONFIND : calmip error flag =',iCalMipErr endif NrMips = nCalMip mipno=min(nCalMip,10) if (nCalMip.GT.10) then write(6,*) 'Warning MUONFIND : > 10 MIPS' endif c Vertex position for track matching : c ------------------------------------ if (COUTAB(VCTVTX).GE.1) then call FETTAB(VCTVTX,ID,1) do i=1,3 vertex(i)=VCTVTX_V(i) enddo else do i=1,3 vertex(i)=0. enddo endif c Set parameters for call to VCTDCA : c ----------------------------------- calerr(1)=20. calerr(2)=20. calerr(3)=20. if (nCalMip.GE.1) then do i=1,min(10,nCalMip) calpos(1)=xCalMip(i) calpos(2)=yCalMip(i) calpos(3)=zCalMip(i) num_match=0 do j=1,200 dist_app(j)=0. vctrhl_index(j)=0. indx(j)=0 enddo c For this MIP, loop over tracks and use VCTDCA : c ----------------------------------------------- if (COUTAB(VCTPAR).GE.1) then do j=1,min(200,COUTAB(VCTPAR)) call FETTAB(VCTPAR,ID,j) vcid=j call vctdca(vcid,calpos,calerr,vctdca_dca, + covprojtrk,covprojcal) c If DCA less than 100 cm, keep this track : c Note that negative DCA is an error condition. c --------------------------------------------- if (vctdca_dca.GE.0..AND.vctdca_dca.LE.100.) then call natrel(vctpar,vctpar_vctrhl,vctrhl,ok) if (ok) then num_match=num_match+1 dist_app(num_match)=vctdca_dca vctrhl_index(num_match)=vctpar_vctrhl c write(6,*) 'vctpar_vctrhl =',vctpar_vctrhl else write(6,*) 'ERROR MUONFIND : VCTPAR w/o VCTRHL entry' endif else if (vctdca_dca.LT.0.) then c write(6,*) 'ERROR MUONFIND : dca =',vctdca_dca endif endif enddo endif c Maximum of 100 matching tracks : c -------------------------------- if (num_match.GT.0) then c Call CERNLLIB routine to sort in ascending order : c -------------------------------------------------- sort_mode=1 sort_nway=0 nsort=0 call SORTZV(dist_app,indx,min(200,num_match),sort_mode, + sort_nway,nsort) c Now, make muon candidate selection cuts : c ----------------------------------------- if (dist_app(indx(1)).LE.dca.AND. + (num_match.EQ.1.OR.dist_app(indx(2)).GE.iso)) then c What is track Pt ? c ------------------ call FETTAB(VCTRHL,ID,vctrhl_index(indx(1))) call cotthetatheta(VCTRHL_tdip,trk_theta) if ((VCTRHL_pgevc*sin(trk_theta)).GE.pt_cut) then c A muon candidate : c ------------------ nomu=nomu+1 c Record Pt for sorting : c ----------------------- pt_sort(nomu)=(VCTRHL_pgevc*sin(trk_theta)) th_sort(nomu)=trk_theta ph_sort(nomu)=VCTRHL_phii endif endif endif enddo endif c Finally, find the two highest momentum muons : c ---------------------------------------------- do i=1,200 indx(i)=0. enddo if (nomu.GE.1) then sort_mode=1 sort_nway=1 nsort=0 call SORTZV(pt_sort,indx,nomu,sort_mode,sort_nway,nsort) do i=1,min(2,nomu) mupt(i)=pt_sort(indx(i)) muth(i)=th_sort(indx(i)) muph(i)=ph_sort(indx(i)) enddo endif return end c