c********************************************************************** c subroutine richdigi c c original design by W. Przygoda (TUM) c last modified on: 01/09/2001 by L. Fabbietti c********************************************************************** implicit none #include "geant321/gcflag.inc" #include "richtups.inc" #include "user.inc" integer nhits, nhmax, nvdim, i, c_hit, p_hit, isect, k integer count(maxmirr) real xsum(maxmirr), ysum(maxmirr) parameter (nhits=11, nhmax=4000, nvdim=3) character*4 iudet integer itrack(nhmax), numvs(nvdim), numbv(nvdim,nhmax) real hits(nhits,nhmax) integer nhret integer cur_tr real x_in, y_in, z_in, theta_in, phi_in, p_in, pxy do k= 1,maxmirr xsum(k) = 0. ysum(k) = 0. count(k) = 0 enddo do isect=1,nsect ! loop over sectors ncer = 0 nhit = 0 nmir = 0 iudet='RDET' ! select detector and sector numvs(1) = 0 numvs(2) = 0 numvs(3) = isect call gfhits('RICH',iudet,nvdim,nhits,nhmax,0,numvs, + itrack,numbv,hits,nhret) IF (NHRET.GT.NHMAX) THEN ! overflow WRITE(6,*) 'ERROR: overflow in richdigi.F (nHits>', & NHMAX,') in RICH',ISECT RETURN ENDIF do i=1,nhret ! loop over hits if(nint(hits(8,i)).eq.50 .and. nint(hits(11,i)).eq.2) then ! a Cerenkov photon reached the pad plane if(ncer.lt.maxckov) then ncer = ncer + 1 xcer(ncer) = lunit*hits(1,i) ! x position in pad plane ycer(ncer) = lunit*hits(2,i) ! y in mm (or cm) ecer(ncer) = hits(7,i) ! photon energy in eV parentcer(ncer) = nint(hits(10,i)) ! track number of parent do k=1,nmircount if(leptrackall(k).eq.parentcer(ncer)) then xsum(k) = xsum(k) + xcer(ncer) ysum(k) = ysum(k) + ycer(ncer) count(k) = count(k) + 1 goto 20 endif enddo 20 continue endif endif if(nint(hits(8,i)).ne.50) then ! direct particle hit if(nint(hits(11,i)).eq.1) then ! particle enters detector cur_tr = nint(hits(10,i)) ! save current track x_in = hits(1,i) ! and entry position (x,y,z) y_in = hits(2,i) z_in = hits(3,i) p_in = hits(7,i)*1.0e3 ! total momentum in MeV/c theta_in = 57.2958*acos(hits(6,i)) ! angle of incidence pxy = sqrt(hits(4,i)**2+hits(5,i)**2) phi_in = 0.0 ! azimuth if(pxy.ne.0.0) phi_in = 57.2958*acos(hits(4,i)/pxy) if(hits(5,i).lt.0.0) phi_in = 360.0 - phi_in else if(cur_tr.eq.nint(hits(10,i))) then ! still same track ? c if(nhit.lt.maxpart .and. hits(9,i).gt.0.0)) then c c ---> take also particles with no energy loss if(nhit.lt.maxpart) then nhit = nhit + 1 xhit(nhit) = lunit*x_in ! x position of incidence yhit(nhit) = lunit*y_in ! y position in mm (or cm) zhit(nhit) = lunit*z_in ! z position parthit(nhit) = hits(8,i) ! particle ID momhit(nhit) = p_in ! total momentum thetahit(nhit) = theta_in ! angle of incidence (0-90 deg) phihit(nhit) = phi_in ! azimuth of incidence (0-360 deg) eloshit(nhit) = hits(9,i) ! energy loss in MeV tlenhit(nhit) = lunit*sqrt((hits(1,i)-x_in)**2 + + (hits(2,i)-y_in)**2 + + (hits(3,i)-z_in)**2) ! path in mm trkhit(nhit) = cur_tr ! GEANT track number else WRITE(6,*) 'ERROR: too many particles (>', & MAXPART,') in RICH',ISECT RETURN endif endif end if enddo ! end loop over hits if(iroot.eq.0) then call hfnt(2000+isect) else call fillrich(isect) endif enddo ! end loop over sectors do k=1,nmircount if (count(k) .gt. 0) then xr(k) = xsum(k) / count(k) yr(k) = ysum(k) / count(k) endif enddo iudet='RMIR' ! select detector and sector numvs(1) = 0 numvs(2) = 0 numvs(3) = 1 call gfhits('RICH',iudet,nvdim,nhits,nhmax,0,numvs, + itrack,numbv,hits,nhret) do i=1,nhret ! loop over hits if(nint(hits(8,i)).eq.2 .or. nint(hits(8,i)).eq.3) then ! e+ or e- if(nint(hits(11,i)).eq.2) then ! particle enters mirror if(nmir.lt.maxpart) then nmir = nmir + 1 xring(nmir) = 0 yring(nmir) = 0 numphot(nmir) = 0 !These 3 variables will be assigned only !if the mirror hit corresponds to a lepton !coming from the radiator; this way leptons !converted in the mirror are not stored. !In fillrich.cc only those leptons are stored !whose number of Cherenkov photons (numphot) !is greater then 0. xlep(nmir) = lunit*nint(hits(1,i)) ! x position of incidence ! on the mirror ylep(nmir) = lunit*nint(hits(2,i)) ! y position in mm (or cm) zlep(nmir) = lunit*nint(hits(3,i)) ! z position lepid(nmir) = nint(hits(8,i)) ! particle ID leptrack(nmir) = nint(hits(10,i)) do k=1,nmircount if (leptrack(nmir) .eq .leptrackall(k)) then xring(nmir) = xr(k) ! average x in pad plane yring(nmir) = yr(k) ! average y numphot(nmir) = numphotall(k) goto 10 endif enddo 10 continue else WRITE(6,*) 'ERROR: too many leptons (>', & MAXPART,') in RICH' RETURN endif endif end if enddo ! end loop over hits call fillrichmirror(1) return end