*********************************************************************** c SUBROUTINE EMCDIGI c c adopted to EMC by K. Lapidus c updated 06/10/2016 R.Holzmann c *********************************************************************** IMPLICIT NONE #include "geant321/gcflag.inc" #include "emctuple.inc" #include "user.inc" #include "kinetups.inc" INTEGER NHDIM, NHMAX, NVDIM PARAMETER (NHDIM=13,NHMAX=18500,NVDIM=2) INTEGER ITRACK(NHMAX), NUMVS(NVDIM), NUMBV(NVDIM,NHMAX) INTEGER NHITS, NTUPID REAL HITS(NHDIM,NHMAX) CHARACTER*4 IUDET CHARACTER CNUM*3, A1*1, A2*2 INTEGER ISECT, I, J, L, NGEATRK REAL NPE_TOT REAL NPE_TRACK(65537) ! 2**16 + 1 REAL SUMPE LOGICAL ISDESCENDANT DO ISECT=1,NSECT ! loop over sectors NTRKS = 0 DO I=1,NMAXCELL ! loop over emc crystals IUDET = 'GLEA' NUMVS(1) = ISECT NUMVS(2) = I CALL GFHITS('EMCG',IUDET,NVDIM,NHDIM,NHMAX,0,NUMVS, & ITRACK,NUMBV,HITS,NHITS) IF (NHITS.GT.NHMAX) THEN ! overflow WRITE(6,*) 'ERROR: overflow in emcdigi.F (nHits>', & NHMAX,') in EMC',ISECT,I RETURN END IF DO J=1,NTRK ! loop over Geant tracks and clear track-wise p.e. counters NPE_TRACK(J) = 0; END DO NPE_TOT = 0 c First loop: integrate the total p.e. yield in a crystal DO L=1,NHITS ! first loop over hits IF (HITS(8,L).GT.0.0) THEN NPE_TRACK(NINT(HITS(7,L))) = & NPE_TRACK(NINT(HITS(7,L))) + HITS(8,L) NPE_TOT = NPE_TOT + HITS(8,L) ENDIF ENDDO ! end first hit loop c make artificial "hit", with total p.e. yield/crystal and crystal nb. in sector IF (NPE_TOT.GT.5.0) THEN ! hard-coded cut of 5 p.e. (~5 MeV) IF (NTRKS.LT.MAXTRKEMC) THEN NTRKS = NTRKS + 1 ! add one output "track/hit" EMCTRK(NTRKS) = -777 ! mark this artificial "track number" EMCDET(NTRKS) = I ! cell number EMCPE(NTRKS) = NPE_TOT ! number of p.e. EMCX(NTRKS) = 0. ! not used EMCY(NTRKS) = 0. ! not used EMCZ(NTRKS) = 0. ! not used EMCTOF(NTRKS) = 0. ! not used EMCMOM(NTRKS) = 0. ! not used EMCLEN(NTRKS) = 0. ! not used ELSE WRITE(6,*) 'ERROR: too many hits (>', & MAXTRKEMC,') in EMC',ISECT,I RETURN END IF END IF c Second loop: Save information relating to incoming Geant tracks. c Take only tracks entering the crystal (from all sides) c i.e. producing a hit within 0.1 mm distance from the borders. DO L=1,NHITS ! second loop over hits c hits(10,l) is INWVOL; hits(12,l) is SAFETY = dist to border IF (HITS(10,L).EQ.1 .AND. HITS(12,L).LT.0.01) THEN c Alternative: Take only tracks bombarding the crystal's front face c IF ( HITS(10,L).EQ.1 .AND. HITS(12,L).LT.0.01 c & .AND. HITS(3,L).LT.-20.99) THEN IF (NTRKS.LT.MAXTRKEMC) THEN NGEATRK = NINT(HITS(7,L)); SUMPE = NPE_TRACK(NGEATRK) DO J=NGEATRK+1,NTRK ! add p.e. of all secondaries of this track IF (ISDESCENDANT(J,NGEATRK)) SUMPE=SUMPE+NPE_TRACK(J) END DO NTRKS = NTRKS + 1 ! add one output "track" EMCTRK(NTRKS) = HITS(7,L) ! Geant track number EMCDET(NTRKS) = HITS(9,L) ! cell number EMCPE(NTRKS) = SUMPE ! p.e. of this Geant track in this module EMCX(NTRKS) = LUNIT*HITS(1,L) ! x pos (in mm) EMCY(NTRKS) = LUNIT*HITS(2,L) ! y pos (in mm) EMCZ(NTRKS) = LUNIT*HITS(3,L) ! z pos (in mm) EMCTOF(NTRKS) = HITS(4,L)*1.0e9 ! mean tof (in ns) EMCMOM(NTRKS) = HITS(13,L)*1.0e3 ! momentum (in MeV/c) EMCLEN(NTRKS) = LUNIT*HITS(11,L) ! track length (in mm) ELSE WRITE(6,*) 'ERROR: too many hits (>', & MAXTRKEMC,') in EMC',ISECT,I RETURN END IF END IF END DO ! end second hit loop END DO ! end emc crystals loop IF (IROOT.EQ.0) THEN CALL HFNT(7000+ISECT) ! 7000 to check ELSE CALL FILLEMC(ISECT) END IF END DO ! end sector loop RETURN END c Check if track ITRK is a descendant of track IGEANTRK c (makes use of track and prevtrack arrays in kinetups.inc) LOGICAL FUNCTION ISDESCENDANT(ITRK, IGEATRK) #include "kinetups.inc" INTEGER ITRK, IGEATRK INTEGER I, IPREV ISDESCENDANT = .FALSE. I = ITRK DO WHILE (.NOT.ISDESCENDANT .AND. I.GT.IGEATRK) c tracks are ordered: desc. have larger track nb than parent IPREV = PREVTRACK(I) IF (IPREV.EQ.IGEATRK) ISDESCENDANT = .TRUE. I = IPREV END DO END