C######################################################################## C C CLEANRHO: function to select clean rho events for calibration purposes C C Author : Heiko Beier C C Selection cuts : - 2 or 3 tracks (dep. on possible e track) C - |eta| < 1.75 & p_t > 150 MeV & opp. charge for tracks C - pt2 < 0.1 (vector sum of e and rho) C - E_3 < 0.2 (highest energy unmatched ConeIsland) C - 0.6 < M_pipi < 1.0 GeV C C all cut values given in data statements below C C######################################################################## C ================================= integer function cleanrho(rhoerr) C ================================= implicit none real etacut(2), ptcut, epzcut, e3cut, pt2cut, mpipicut(2) data etacut /1.75, -1.75/ data ptcut /0.15/ data epzcut /40.0/ data e3cut /0.2/ data pt2cut /0.1/ data mpipicut /0.6,1.0/ integer rhoerr, ierr integer rho_trkcuts, rho_bkgcuts cleanrho = 0 rhoerr = 0 * DIS event? BPC or SRTD/CAL? call rho_elec(ierr) if (ierr.ne.0) then rhoerr = 1 ! no (unique) electron identified return endif * variables related to matching call zes_matsch(ierr) if (ierr.ne.0) then rhoerr = 2 ! some error with matching return endif * check tracking variables if (rho_trkcuts(etacut,ptcut,epzcut).eq.0) then rhoerr = 3 ! failed tracking cuts return endif * apply background cuts if (rho_bkgcuts(e3cut,pt2cut,mpipicut).eq.0) then rhoerr = 4 ! failed background cuts return endif cleanrho = 1 return end C ========================= subroutine rho_elec(ierr) C ========================= implicit none #include "zescommon.inc" #include "zestmpcommon.inc" #include "rhowork.inc" integer ierr, cflag, sflag, bflag REAL xv(3), tmp_x, tmp_y, tmp_z ierr = 0 if (ZES_Eeu_si.gt.0.) then call calbox(cflag) endif if (ZES_Esrtd_si.gt.0.) then call srtdbox(sflag) endif if (ZES_Eebpc.gt.0.) then call bpcbox(bflag) endif if (bflag.eq.0.and.cflag.eq.0.and.sflag.eq.0) then ierr = 1 return endif if (bflag.eq.1.and.(sflag.eq.1.or.cflag.eq.1)) then C print *,'Electron candidate in both BPC and SRTD/CAL!' ierr = 2 return endif if (sflag.eq.1) then CALL SETVTX(XV) ebest = ZES_Eeu_si tmp_x = ZES_xsrtd_si - XV(1) tmp_y = ZES_ysrtd_si - XV(2) tmp_z = -147.9 - XV(3) tbest = ATan2(sqrt(tmp_x**2 + tmp_y**2), & tmp_z) pbest = atan2(tmp_y,tmp_x) elseif (cflag.eq.1) then ebest = ZES_ecorr_si tbest = ZES_theta_si pbest = ZES_phi_si elseif (bflag.eq.1) then ebest = ZES_Eebpc tmp_x = ZES_pxbpc tmp_y = ZES_pybpc tmp_z = ZES_pzbpc tbest = ATan2(sqrt(tmp_x**2 + tmp_y**2), & tmp_z) pbest = atan2(tmp_y,tmp_x) endif return end C ========================= subroutine calbox(cflag) C ========================= implicit none #include "zescommon.inc" INTEGER CFLAG REAL xcut, ycut DATA xcut, ycut /16.0,12.0/ IF (ABS(ZES_Xeu_si).GT.xcut .OR. & ABS(ZES_Yeu_si).GT.ycut) THEN cflag = 1 !!! accept ELSE cflag = 0 !!! reject ENDIF RETURN END C ========================== subroutine srtdbox(sflag) C ========================== implicit none #include "zescommon.inc" INTEGER sFLAG REAL xcut, ycut DATA xcut, ycut /13.0,8.0/ IF (ABS(ZES_Xsrtd_si).GT.xcut .OR. & ABS(ZES_Ysrtd_si).GT.ycut) THEN sflag = 1 !!! accept ELSE sflag = 0 !!! reject ENDIF RETURN END C ======================== subroutine bpcbox(bflag) C ======================== C 95 fiducial cut definition from Teresa. C implicit none #include "zescommon.inc" INTEGER bFLAG IF (ZES_xebpc.LE. 5.17 .OR. & ZES_xebpc.GE. 9.97 .OR. & ZES_yebpc.LE.-2.82 .OR. & ZES_yebpc.GE. 3.18) THEN bflag = 0 !!! reject ELSE bflag = 1 !!! accept ENDIF RETURN END C =========================== subroutine zes_matsch(ierr) * 29/5/98 A. Bertolin changed float(PROBCUT) in PROBCUT twice C =========================== implicit none #include "rhowork.inc" #include "vmvars.inc" #include "zisles.inc" ! /* ! Island Common from ZUFOs, gives all clustering info */ INTEGER IsItEcone real PROBCUT, DMINC Parameter (PROBCUT = 10.) C Parameter (DMINC = 0.) integer NISLMAX Parameter (NISLMAX=50) real maxprob, matprob, dist, disterr, mat(NISLMAX), & mprob(NISLMAX) real islpos(3,NISLMAX), islth(NISLMAX) integer materr, ierr, icon, itrk, first data first /1/ DMINC = 0. e3 = 0. ierr = 0 if (first) then first = 0 if (IslTyp(1).ne.11) then ! IslTyp variable from common ZIsles. WRITE(*,*)'Warning zes_matsch: ConeIsls reconstructed with' WRITE(*,*)' IslTyp =',IslTyp(1),zcIsl_Mode WRITE(*,*)' MATSCH expects IslTyp = 11' endif endif if (nIsl.gt.NISLMAX) then ierr = 1 C print * C & ,'# ConeIslands exceeds maximum specified in zes_matsch()!' return endif do icon=1,NIsl ! NIsl from ZES_ZUFOS common maxprob = 0. islpos(1,icon) = xIsl(icon) islpos(2,icon) = yIsl(icon) islpos(3,icon) = zIsl(icon) islth(icon) = atan2(sqrt(xIsl(icon)**2.+yIsl(icon)**2.), & zIsl(icon)) * inner loop over all tracks do itrk=1,ntrvtx if (IsItEcone(icon).eq.1) then * Matsch V2.03 for hadrons call matsch(itrk,1,islpos(1,icon),islth(icon), & eIsl(icon),emcEIsl(icon)/eIsl(icon),1,0, & DMINC,dist,disterr,matprob,materr) else * Matsch V2.03 for electrons call matsch(itrk,1,islpos(1,icon),islth(icon), & eIsl(icon),emcEIsl(icon)/eIsl(icon),2,0, & DMINC,dist,disterr,matprob,materr) endif if (materr.eq.0) then if (matprob.gt.maxprob) then maxprob = matprob mat(icon) = itrk mprob(icon) = maxprob endif endif enddo ! loop over tracks if (mat(icon).gt.0) then * with MatschV2 apply prob cut to both pions and elec clusters if (IsItEcone(icon).eq.1) then if (mprob(icon).gt.(PROBCUT/1000.)) then ideltrk = mat(icon) ! mark electron track endif endif * get E3 for hadronic clusters if (mprob(icon).lt.(PROBCUT/1000.).and. & mat(icon).ne.ideltrk ) then if (eIsl(icon).gt.e3) then e3=eIsl(icon) endif endif endif enddo ! loop over clusters return end C ================================================= integer function rho_trkcuts(etacut,ptcut,epzcut) C ================================================= implicit none #include "rhowork.inc" #include "vmvars.inc" ! /* ! vector meson specific info from Teresa's VMSEL routine */ real etacut(2), ptcut, epzcut real yjb, ee Parameter (ee=27.5) rho_trkcuts = 0 * # vertex tracks if (ntrvtx.lt.2.or.ntrvtx.gt.3) return * # tracks compatible with electron track? if (ntrvtx.eq.2.and.ideltrk.ne.0) return if (ntrvtx.eq.3.and.ideltrk.eq.0) return * cuts on eta, pt (pion tracks only) if ( eta(1).gt.etacut(1).or. & eta(1).lt.etacut(2).or. & eta(2).gt.etacut(1).or. & eta(2).lt.etacut(2).or. & sqrt(ptr(1,1)**2.+ptr(1,2)**2.).lt.ptcut.or. & sqrt(ptr(2,1)**2.+ptr(2,2)**2.).lt.ptcut ) then return endif * cut on E-pz (calculated using tracks) yjb = 1./2./ee * (ptr4pi(1)-ptr(1,3) & + ptr4pi(2)-ptr(2,3) ) if ( (2.*ee*yjb + ebest*(1.-cos(tbest))).lt.epzcut) return * opposite sign if (charge(1)*charge(2).eq.1) return * passed rho_trkcuts = 1 return end C =================================================== integer function rho_bkgcuts(e3cut,pt2cut,mpipicut) C =================================================== implicit none #include "rhowork.inc" #include "vmvars.inc" ! /* ! vector meson specific info from Teresa's VMSEL routine */ #include "zescommon.inc" real e3cut, pt2cut, mpipicut(2) real ee, eecs, pt2 Parameter (ee=27.5) rho_bkgcuts = 0 * cut on unmatched energy if (e3.gt.e3cut) return * cut on pt2 eecs = (2.*ee - ((ptr4pi(1)+ptr4pi(2))-(ptr(1,3)+ptr(2,3)))) / & (1. - cos(tbest)) pt2 = ((ptr(1,1)+ptr(2,1)) + eecs*sin(tbest)*cos(pbest))**2. + & ((ptr(1,2)+ptr(2,2)) + eecs*sin(tbest)*cos(pbest))**2. if (pt2.gt.pt2cut) return * cut on mass if (ZES_Mtrknoe_pi.lt.mpipicut(1).or. & ZES_Mtrknoe_pi.gt.mpipicut(2)) return * passed rho_bkgcuts = 1 return end C ================================= integer function IsItEcone(Icone) C ================================= IMPLICIT NONE #include "zescommon.inc" #include "zestmpcommon.inc" #include "zisles.inc" INTEGER Iflag, I, J, Icone Iflag = 0 DO I=1,ZESTMP_Ncell_si DO J=1,NrcIsl(Icone) IF (ZESTMP_CellList_si(I).EQ.pCellIsl(J,Icone)) Iflag=1 ENDDO ENDDO IsItEcone = Iflag IF (Iflag.eq.1) WRITE(*,*) 'electron found' RETURN END