* * $Id: gustep.F,v 1.1.1.1 1995/10/24 10:21:46 cernlib Exp $ * * $Log: gustep.F,v $ * Revision 1.1.1.1 1995/10/24 10:21:46 cernlib * Geant * * #include "geant321/pilot.h" *CMZ : 3.21/02 29/03/94 15.41.25 by S.Giani *-- Author : SUBROUTINE GUSTEP C. C. ****************************************************************** C. * * C. * User routine called at the end of each tracking step * C. * INWVOL is different from 0 when the track has reached * C. * a volume boundary * C. * ISTOP is different from 0 if the track has stopped * C. * * C. * ==>Called by : GTRACK * C. * * C. ****************************************************************** C. #include "geant321/gcflag.inc" #include "geant321/gckine.inc" #include "geant321/gcking.inc" #include "geant321/gctmed.inc" #include "geant321/gctrak.inc" #include "geant321/gcvolu.inc" C. C. ------------------------------------------------------------------ C. real cost, sint, phi real cosp, sinp character*4 cnames(15),cnamec(15),cnamec1(15) equivalence (names(1),cnames(1)) equivalence (namec(1),cnamec(1)) equivalence (namec1(1),cnamec1(1)) integer ncaptured,ndetected common /fiberstats/ncaptured,ndetected data ncaptured,ndetected/0,0/ real rndm(10) real x0(6), x1(3) integer medi0, medi1 real u(3), pdotu integer ierr integer identry,idloss common /saveme/identry,idloss data identry,idloss/0,0/ real rindex(2) real rabsco real cladthick real klong,kappa real lambda data lambda/485e-7/ c If it is a geantino, stop it and convert it to a Cerenkov c photon of a random direction and energy. This is because c Geant3 only supports injection of optical photons into the c simulation through secondary interactions. if (ipart .eq. 48) then ngphot = 1 xphot(1,1) = vert(1) xphot(2,1) = vert(2) xphot(3,1) = vert(3) call grndm(rndm,2) cost = rndm(1)*2-1 sint = sqrt(1-cost**2) phi = rndm(2)*2*3.1415926 xphot(4,1) = sint*cos(phi) xphot(5,1) = sint*sin(phi) xphot(6,1) = cost xphot(7,1) = pvert(4) xphot(8,1) = cost*cos(phi) xphot(9,1) = cost*sin(phi) xphot(10,1) = -sint xphot(11,1) = tofg call gskpho(0) istop = 1 endif c If it is an optical photon, and it is reflecting from an interface c between two dielectrics, take into account the absorption of the c evanescent wave in the outer cladding, and don't just assume that c total internal reflection is impervious to the conditions of the c outer cladding layers. if (ipart.eq.50 .and. nmec.gt.1 .and. lmec(nmec).eq.106) then x0(1) = vect(1) x0(2) = vect(2) x0(3) = vect(3) x0(4) = vect(4) x0(5) = vect(5) x0(6) = vect(6) medi0 = numed if (medi0 .eq. 100) then rindex(1) = 1.60 rindex(2) = 1.49 cladthick = 0.008 rabsco = 1e5 medi1 = 101 elseif (medi0 .eq. 101) then rindex(1) = 1.49 rindex(2) = 1.42 cladthick = 0.004 rabsco = 1e5 medi1 = 102 else rindex(1) = 1.42 rindex(2) = 1.0 cladthick = 0.0001 rabsco = 1 medi1 = 1 endif x1(1) = vect(1)+vect(4)*0.01 x1(2) = vect(2)+vect(5)*0.01 x1(3) = vect(3)+vect(6)*0.01 call glisur(x0, x1, medi0, medi1, u, pdotu, ierr) if (ierr .eq. 0) then klong = sqrt(1-pdotu**2) if (klong .gt. rindex(2)/rindex(1)) then kappa = sqrt(klong**2 - (rindex(2)/rindex(1))**2) absprob = rabsco*exp(-(2*3.1416/lambda)*kappa*cladthick) call grndm(rndm,1) if (rndm(1) .lt. absprob) then * print *, 'absorption',absprob,' in ',cnames(nlevel) istop = 1 else * print *, 'reflection',absprob,' in ',cnames(nlevel), * + kappa,(2*3.1416/lambda)*kappa*cladthick continue endif endif endif endif c Do the detection of trapped and transmitted photons, c correcting for those that are falsely reflected from the c unphysical surfaces between connected arc and straight sections. if (nlevel .gt. 1) then if (inwvol .eq. 1 .and. cnames(2) .eq. 'OCA1' .and. + vect(4) .gt. 0 .and. abs(vect(1) - 10) .lt. 0.1) then if (identry .ne. idevt) then ncaptured = ncaptured + 1 endif identry = idevt else if (inwvol .eq. 1 .and. cnames(2) .eq. 'OCS1' .and. + vect(4) .lt. 0 .and. abs(vect(1) - 10) .lt. 0.1) then if (identry .eq. idevt .and. idloss .ne. idevt) then ncaptured = ncaptured - 1 endif idloss = idevt identry = 0 else if (cnames(2) .eq. 'OCS2' .and. vect(1) .gt. 170.) then if (cnames(nlevel) .eq. 'COS2') then print *, 'cos(theta) =', vect(4) ndetected = ndetected + 1 endif istop = 1 endif endif #if DEBUGGING_STUFF if (inwvol .eq. 2 .and. cnames(nlevel) .eq. 'COA1') then x0(1) = vect(1) x0(2) = vect(2) x0(3) = vect(3) x0(4) = vect(4) x0(5) = vect(5) x0(6) = vect(6) medi0 = numed else if (inwvol .eq. 1 .and. cnames(nlevel) .eq. 'COA2' .and. + medi0 .gt. 0) then x1(1) = vect(1) x1(2) = vect(2) x1(3) = vect(3) medi1 = numed call glisur(x0, x1, medi0, medi1, u, pdotu, ierr) print *, 'transit from medium',medi0,' to',medi1,' at',x1 print *, ' unit normal to surface is',u print *, ' new direction is',vect(4),vect(5),vect(6) else if (medi0 .ne. 0) then print *, 'leaving medium',medi0,' for medium',numed medi0 = 0 endif if (abs(vect(3)) .lt. 2.0 .and. + abs(vect(2)+2.33) .lt. 2.0 .and. + abs(vect(1)-16.45) .lt. 2.0) then print *, 'new direction:',vect(4),vect(5),vect(6), + ' in ',cnames(nlevel) print *, ' and epsil is',epsil endif #endif call gdebug END