* #define _MAXEVE_ 1000 * #define _MAXTOTHITS_ 35*_MAXEVE_ #if _MAXEVE_ > 100000 #error "Number of event to read from library exceeds library size" #endif * subroutine addhits_lib(x0,y0,nx,ny,e,chopt) implicit none #include "cphoto.inc" #include "adcgam_bk.inc" * character*(1) chopt real x0, y0, nx, ny, e * integer neve parameter(neve = _MAXEVE_) * integer ech common/ech_common/ ech(0:1,-MCOL:MCOL,-MROW:MROW) * integer MAXEVE, MAXTOTHITS, NLIBEXT parameter(MAXEVE = _MAXEVE_, MAXTOTHITS = _MAXTOTHITS_) parameter(NLIBEXT = 105) * integer stored_eve_addr integer*1 stored_eve_len integer*2 stored_hit_val common /lib_store_com/ & stored_eve_addr(NLIBEXT,MAXEVE), & stored_eve_len(NLIBEXT,MAXEVE), & stored_hit_val(NLIBEXT,MAXTOTHITS) * integer hit_ind, i, j, iaddr, iid, ie, blklen, inum, ien, ihole, & kx, ky, kx1, kx2, ky1, ky2, iquadrant, ixyind, & ix, iy, nhit integer iseq(NLIBEXT) save iseq real x, y, dx, dy, dz, rndm, shower_depth, ecell logical first data first /.true./ save first * integer idim parameter(idim=3) integer a(-idim:idim,-idim:idim) * integer*2 i2word character*(128) path, fnam * integer nebin parameter(nebin = 7) real table_energy(nebin) data table_energy /0.1, 0.25, 0.5, 1., 2., 4., 8./ save table_energy * real eratio, esum, escals, f, eshower_mean(nebin) data eshower_mean /.902, .928, .936, .941, .953, .953, .953/ save eshower_mean * if (first) then first = .false. call getenv('LIB_DIR',path) * if (chopt.eq.'r') then ! pick up random event print*, 'start reading library randomly' elseif (chopt.eq.'s') then ! pick up sequential event print*, 'start reading library sequentially' else print*, 'read_event_from_lib error: unknwn chopt value ', & chopt call exit(1) endif * write(*,'(a1,$)') '|' do inum = 1, NLIBEXT iseq(inum) = 0 * fnam = 'lib.dat.' write(fnam(9:12),'(I4.4)') inum write(*,'(a1,$)') '*' call read_lib_file(inum, & path(:blklen(path))//'/'//fnam(:blklen(fnam)), neve) enddo write(*,'(a1)') '|' * endif * if(igeometry.ne.1.or.icase.ne.1) then ! if not hybrid print*, 'Shower library is not implemented for LG yet' call exit(1) endif * i_gen_method = 2 ! read from library dz = shower_depth(icase,e) x = x0 + nx*dz y = y0 + ny*dz kx = nint(x/(dtarg(icase,1)*2.)) ky = nint(y/(dtarg(icase,2)*2.)) * if(icase.eq.1) then ! PWO if(.not.(abs(x0).lt.60. .and. abs(y0).lt.60.)) return ihole = 3 if(iabs(kx).gt.NCOL(1).or.iabs(ky).gt.NROW(1)) return ! out of calorimeter if(iabs(kx).lt.ihole .and. iabs(ky).lt.ihole) return ! hole else if(icase.eq.0) then ! LG if(abs(x0).lt.60. .and. abs(y0).lt.60.) return ihole = 2 if(kx*kx+ky*ky.gt.960) return ! out of calorimeter if(iabs(kx).lt.ihole .and. iabs(ky).lt.ihole) return ! hole else ihole = 0 endif * call get_lib_energy(ien,e,escale) * dx = x/(dtarg(icase,1)*2.)-kx dy = y/(dtarg(icase,2)*2.)-ky * call get_cell_sector(iquadrant,ixyind,dx,dy,xmislib,ymislib) * inum = (ien-1)*15+ixyind iseq(inum) = iseq(inum) + 1 * c print*, 'ie ixy iqua dx dy', ien, ixyind, iquadrant, dx, dy * if (chopt.eq.'r') then ! pick up random event hit_ind = 1 + int(rndm(-1.)*neve) elseif (chopt.eq.'s') then ! pick up sequential event hit_ind = 1 + mod(iseq(inum)-1,neve) else print*, 'read_event_from_lib error: unknwn chopt value ', chopt call exit(1) endif * iaddr = stored_eve_addr(inum,hit_ind) * nhit = stored_eve_len(inum,hit_ind) call vzero(a,(2*idim+1)**2) esum = 0. do i = 1, nhit i2word = stored_hit_val(inum,iaddr+i) iid = iand(i2word,63) ie = ishft(i2word,-6) * if(iid.eq.0) then ix = 0 iy = 0 else if(iid.gt.24) iid = iid + 1 ix = idim - mod((iid-1),(2*idim+1)) iy = idim - (iid-1)/(2*idim+1) endif esum = esum + ie enddo * eratio = 0.01*esum/table_energy(ien) if(ien.eq.7) then f = 0.70/sqrt(sqrt(escale)) else f = 0.75/sqrt(sqrt(escale)) endif escals = f + eshower_mean(ien)*(1.-f)/eratio * do i = 1, nhit i2word = stored_hit_val(inum,iaddr+i) iid = iand(i2word,63) ie = ishft(i2word,-6) * if(iid.eq.0) then ix = 0 iy = 0 else if(iid.gt.24) iid = iid + 1 ix = idim - mod((iid-1),(2*idim+1)) iy = idim - (iid-1)/(2*idim+1) endif ecell = ie*escale*escals ecell = ecell * 1.04853 * (1. + 0.03*(escale-1.)) * a(ix,iy) = nint(100*ecell) enddo * call rotate_shower(iquadrant,a) * kx1 = max(-NCOL(icase),kx-idim) kx2 = min( NCOL(icase),kx+idim) ky1 = max(-NROW(icase),ky-idim) ky2 = min( NROW(icase),ky+idim) * do i = kx1, kx2 do j = ky1, ky2 if(icase.eq.0.and.i*i+j*j.ge.961) goto 120 ! out of LG calorimeter if(iabs(i).lt.ihole .and. iabs(j).lt.ihole) goto 120 ! hole ech(icase,i,j) = ech(icase,i,j) + a(i-kx,j-ky) 120 continue enddo enddo * * return end * * -------------------------------------- * subroutine read_lib_file(inum,fnam,neve) implicit none * integer inum, neve character*(64) fnam * integer mhit, liblun parameter(mhit = 64, liblun = 21) * integer pos, i, irec, ie, nint, ihit, ihit_addr, ihit_addr_old integer*1 rec_len, id_list(mhit), i1(4) integer*2 hitw ! hit word: mixtured id and energy real fe_list(mhit), few equivalence (few,i1) * integer MAXEVE, MAXTOTHITS, NLIBEXT parameter(MAXEVE = _MAXEVE_, MAXTOTHITS = _MAXTOTHITS_) parameter(NLIBEXT = 105) * integer stored_eve_addr integer*1 stored_eve_len integer*2 stored_hit_val common /lib_store_com/ & stored_eve_addr(NLIBEXT,MAXEVE), & stored_eve_len(NLIBEXT,MAXEVE), & stored_hit_val(NLIBEXT,MAXTOTHITS) * if(neve.gt.MAXEVE .or. inum.lt.1 .or. inum.gt.NLIBEXT) then print*, 'read_lib_file error: bad values for inum or neve ', & inum, neve call exit(1) endif * pos = 1 open(liblun, file = fnam, form = 'unformatted', & access = 'direct', status = 'old', recl = 1) * ihit = 0 ihit_addr = 0 do irec = 1, neve read(liblun, REC = pos) rec_len pos = pos + 1 ihit = ihit + 1 * if(ihit_addr+rec_len.gt.MAXTOTHITS) then print*, 'read_lib_file error: too many hits in library' call exit(1) endif * ihit_addr_old = ihit_addr stored_eve_addr(inum,ihit) = ihit_addr * do i = 1, rec_len read(liblun, REC = pos) id_list(i) pos = pos + 1 read(liblun, REC = pos) i1(1) pos = pos + 1 read(liblun, REC = pos) i1(2) pos = pos + 1 read(liblun, REC = pos) i1(3) pos = pos + 1 read(liblun, REC = pos) i1(4) pos = pos + 1 * fe_list(i) = few * ie = nint(fe_list(i)*100.) if(ie.eq.0) goto 50 ! eliminate chan. if below threshold * if(id_list(i).lt.0 .or. id_list(i).gt.60 .or. & ie.lt.0 .or. ie.gt.1000) then print*, 'read_lib_file error: bad hit values ', & id_list(i), fe_list(i) call exit(1) endif hitw = iand(id_list(i),63) + ishft(ie,6) ihit_addr = ihit_addr + 1 stored_hit_val(inum,ihit_addr) = hitw 50 continue enddo stored_eve_len(inum,ihit) = ihit_addr - ihit_addr_old enddo close(liblun) cc print*, 'last read record address = ', ihit_addr * return end * * -------------------------------------- * subroutine get_lib_energy(ien,e,escale) implicit none integer ien real e, escale * integer nebin parameter(nebin = 7) real table_energy(nebin) data table_energy /0.1, 0.25, 0.5, 1., 2., 4., 8./ save table_energy integer i real ecut * do i = 1, nebin-1 ecut = sqrt(table_energy(i)*table_energy(i+1)) if(e.lt.ecut) then ien = i goto 110 endif enddo ien = nebin 110 continue escale = e/table_energy(ien) * return end * * -------------------------------------- * subroutine get_cell_sector(iq,it,x,y,xmis,ymis) implicit none * integer iq, it ! quadrant number and hit position index real x, y ! in cell size units real xmis, ymis ! misdistances due to lib. digitization integer ldim parameter(ldim = 4) integer ix, iy, i, nint, & is(-ldim:ldim,-ldim:ldim), & ie(0:ldim,0:ldim) data is / 5, 6, 6, 6, 6, 7, 7, 7, 7, & 5, 5, 6, 6, 6, 7, 7, 7, 8, & 5, 5, 5, 6, 6, 7, 7, 8, 8, & 5, 5, 5, 5, 6, 7, 8, 8, 8, & 4, 4, 4, 4, 1, 1, 1, 1, 1, & 4, 4, 4, 3, 2, 1, 1, 1, 1, & 4, 4, 3, 3, 2, 2, 1, 1, 1, & 4, 3, 3, 3, 2, 2, 2, 1, 1, & 3, 3, 3, 3, 2, 2, 2, 2, 1/ data ie / 1, 2, 4, 7, 11, & 0, 3, 5, 8, 12, & 0, 0, 6, 9, 13, & 0, 0, 0,10, 14, & 0, 0, 0, 0, 15/ save is, ie * ix = nint(2*ldim*x) iy = nint(2*ldim*y) xmis = x - ix/float(2*ldim) ymis = y - iy/float(2*ldim) * if(iabs(ix).le.ldim.and.iabs(iy).le.ldim) then iq = is(ix,iy) else print*, 'Error: array out of range', x, y iq = 0 call exit(1) endif * ix = iabs(ix) iy = iabs(iy) if(iy.gt.ix) then i = ix ix = iy iy = i endif it = ie(ix,iy) * return end * * -------------------------------------- * subroutine rotate_shower(ics,a) implicit none integer ics integer idim parameter(idim=3) integer i, j integer a(-idim:idim,-idim:idim) integer buf logical trans, reflectx, reflecty * trans = ics.eq.2 .or. ics.eq.3 .or. ics.eq.6 .or. ics.eq.7 reflectx = ics.eq.3 .or. ics.eq.4 .or. ics.eq.5 .or. ics.eq.6 reflecty = ics.eq.5 .or. ics.eq.6 .or. ics.eq.7 .or. ics.eq.8 * if(trans) then do i = -idim, idim-1 do j = i+1, idim buf = a(i,j) a(i,j) = a(j,i) a(j,i) = buf enddo enddo endif * if(reflectx) then do i = -idim, -1 do j = -idim, idim buf = a(i,j) a(i,j) = a(-i,j) a(-i,j) = buf enddo enddo endif * if(reflecty) then do i = -idim, idim do j = -idim, -1 buf = a(i,j) a(i,j) = a(i,-j) a(i,-j) = buf enddo enddo endif * return end * * -------------------------------------- * Integer Function BLKLEN(Text) Implicit none Character*(*) Text Integer I,J,K I = len(Text) Do J = 1,I K = J If (ichar(text(j:j)).eq.0) Goto 1 EndDo Goto 2 1 Continue I = K - 1 2 Continue Do J = I,1,-1 K = J If (text(j:j).ne.' ') Goto 3 EndDo 3 Continue If (K.eq.1 .and. (text(1:1).eq.' ' .or. & ichar(text(1:1)).eq.0)) K = 0 blklen = K Return End *