*********************************************************************** c SUBROUTINE EMCSTEP c c adapted to EMC by K.Lapidus c updated 05/10/2016 R.Holzmann c *********************************************************************** IMPLICIT NONE #include "geantdef.h" #include "geant321/gcflag.inc" #include "geant321/gckine.inc" #include "geant321/gcking.inc" #include "geant321/gclist.inc" #include "geant321/gcsets.inc" #include "geant321/gctmed.inc" #include "geant321/gctrak.inc" #include "geant321/gcvolu.inc" #include "emctuple.inc" #include "geamedia.inc" INTEGER NH PARAMETER (NH=13) REAL HITS(NH) REAL XM(3), XD(3), XDC(3), PM(3), PD(3) REAL XU(3) REAL PMTEFF, RNDM, EFFI, DIVDIF CHARACTER*4 CHSET EQUIVALENCE (IHSET,CHSET) INTEGER I, J, IHIT INTEGER DETNR, CELLINDEX, PMTTYPE INTEGER*8 HIST1, HIST2 REAL NPE ! number of photoelectrons (p.e.) SAVE NPE REAL RES IF(CHSET.NE.'EMCG') RETURN ! particle is not in the EMC IF(IPART.EQ.4) RETURN ! skip neutrinos c DETNR = IDTYPE ! not used c accumulate photo electrons produced in the current step IF(NUMED.EQ.NUMED_LEADGLASS) THEN ! lead glass IF(INWVOL.EQ.1) NPE = 0.0 ! track entered, or a new track IF(NGPHOT.GT.0) THEN c retrieve PMT type and set pointer to corresonding histograms CELLINDEX = (NUMBER(3) - 1) * NMAXCELL + NUMBER(4) PMTTYPE = PMTLOOKUP(CELLINDEX) IF (PMTTYPE.EQ.2) THEN HIST1 = POINT3 HIST2 = POINT4 ELSE HIST1 = POINT1 HIST2 = POINT2 ENDIF DO I=1,NGPHOT ! loop over generated Cherenkov photons PMTEFF = RNDM(1) c linear interpolation in eff table using CERNLIB function divdif (E105) EFFI = DIVDIF(EMCPHOTEFF,EMCPHOTLEN,EMCOPTBIN, & 1.E9*XPHOT(7,I),1) IF (PMTEFF.LT.EFFI) THEN ! this Cherenkov photon would be converted into p.e. DO J=1,3 XM(J) = XPHOT(J,I) PM(J) = XPHOT(J+3,I) ENDDO CALL GMTOD(XM,XD,1) ! transform coordinates into detector frame CALL GMTOD(PM,PD,2) RES = 0. c call external code PROFAST to look up propagation probability of Cherenkov photon CALL PROFAST(HIST1,HIST2,XD(1),XD(2),XD(3), & PD(1),PD(2),PD(3),1.E9*XPHOT(7,I),RES) IF (RES.GT.0.0) NPE = NPE + RES ! add this amount of p.e. ENDIF ENDDO ! end loop over generated Cerenkov photons ENDIF ENDIF IF(INWVOL.EQ.2 .OR. INWVOL.EQ.1 .OR. ISTOP.GE.1) THEN c ! volume left or entered, or particle stopped or decayed. DO I=1,3 XM(I) = VECT(I) ENDDO CALL GMTOD(XM,XD,1) HITS(1) = XD(1) ! x of hit in detector coordinates HITS(2) = XD(2) ! y " HITS(3) = XD(3) ! z " HITS(4) = TOFG ! tof of current hit HITS(5) = FLOAT(IPART) ! GEANT particle ID HITS(6) = NUMBER(3) ! volume copy number at level 3 (=sector) HITS(7) = FLOAT(ITRA) ! track nb. HITS(8) = NPE ! p.e. yield from this track HITS(9) = NUMBER(4) ! volume copy number at level 4 (=cell) HITS(10) = FLOAT(INWVOL) ! tracking flag (in common GCTRAK) HITS(11) = SLENG ! length of current track HITS(12) = SAFETY ! distance to the border HITS(13) = VECT(7) ! total momentum CALL GSAHIT(ISET,IDET,ITRA,NUMBV,HITS,IHIT) ! store hit for emcdigi ENDIF RETURN END