PROGRAM Main C This example makes the ray-tracing procedure for a single C user defined ray within the user defined Winston Cone IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X0(2),V0(3),X1(2),V1(3) REAL*4 RANO REAL SINT,COST,SINP,COSP,THETA(1000000),PHI(1000000),th REAL pi180,xps(1000000),yps(1000000) REAL px,py,pz,px1,py1,pz1,th1,ph1 C REAL px(1000000),py(1000000),pz(1000000) C integer mhbook,NOENT C parameter (mhbook=100000) C C make index integer C integer c integer nsuccess c c Updated program to get input from the 90deg light guides for the c outer layers of the Bcal. Elton Smith 2/12/08 c c assume input histogram covers a circle with a radius of 2.8284 cm. c TIN=2*2.8284 ! entrance diameter of the Winston Cone c TOUT=2*2.54 ! exit diameter of the Winston Cone c c TIN=2*2.625 ! entrance diameter of the Winston Cone outer c TOUT=2*1.3 ! exit diameter of the Winston Cone outer c c TIN=2*1.685 ! entrance diameter of the Winston Cone inner fine mesh TIN=2*1.35 ! entrance diameter of the Winston Cone inner fine mesh necked down to 1.35 TOUT=2*0.875 ! exit diameter of the Winston Cone inner fine mesh c c TIN=2*1.0 ! entrance diameter of the Winston Cone inner SiPM c TOUT=2*0.6 ! exit diameter of the Winston Cone inner RR=1.49 ! ratio between internal and external C refractive indeces C RR=0. ! if RR.le.1 100% reflection independent of C incidence angle NMAX=100 ! maximum allowed # of bounces for a single ray nsuccess = 0 ! initialize number of photons that pass CALL WICOIN(TIN,TOUT,RR,NMAX,ELLE) PRINT*,' Winston Cone LENGTH =',ELLE C Initialize Histograms * call GZEBRA(mzebra) * call HLIMIT(mhbook) CALL HBOOKINIT CALL HBOOKSTART call HROPEN(2,'BCAL','gntbcal.hbook',' ',4096,istat) print *,' Input Hbook file is gntbcal.hbook' print *,' Input Radius =',TIN/2 print *,' Output Radius =',TOUT/2 print *,' Index of Refraction=',RR CALL HIDOPT(0,'STAT') if (istat.ne.0) then stop 'cry foul' endif c c numbering equivalent for Elton's geant output: c hist old id elton id c lg CPC face 102 214 c exit phi 302 304 c exit theta 305 305 call HRIN(214,999,0) call HNOENT(214,NEVENT) print *,' Number of events in histogram 214 used=',nevent * call HRIN(105,999,0) * call HRIN(106,999,0) * call HRIN(107,999,0) call HRIN(304,999,0) call HRIN(305,999,0) * noent=10000 c c for straight WC, eliminate offset z0out = 0. c z0out = 7.125 do i=1,NEVENT c c call random variables for light ray angle and position and store in array c xps(i)=0. yps(i)=0. 100 CALL HRNDM2(214,xps(i),yps(i)) c c translate to 0, 0 c xps(i) = xps(i) - z0out c print *,' event=',i,' X0(1)=',xps(i),' X0(2)=',yps(i) c c check for generation outside radius due to histogram granularity c if (SQRT(xps(i)**2+yps(i)**2).ge.TIN/2) GOTO 100 102 THETA(i)=HRNDM1(305) if(theta(i).gt.90) goto 102 PHI(i)=HRNDM1(304) * px(i)=HRNDM1(106) * py(i)=HRNDM1(107) * 101 pz(i)=HRNDM1(105) * if(pz(i).lt.0) GOTO 101 c if ((i/10)*10 .eq. i) then c print *, ' Event=',i,' phi=',phi(i),' theta=',theta(i) c print *, ' Event=',i,' xps=',xps(i),' yps=',yps(i) endif enddo c call HDELET(214) * call HDELET(105) * call HDELET(106) * call HDELET(107) call HDELET(304) call HDELET(305) call HREND('BCAL') close(2) C------- loop on input data * C=0 * CMAX=10 * DO WHILE ( C .LT. CMAX ) call HBOOK2 (3,'Input y vs x',100,-4.,4.,100,-4.,4.,0.) CALL HBOOK1 (4,'Input px',100,-1.,1.,0.) CALL HBOOK1 (5,'Input py ',100,-1.,1.,0.) CALL HBOOK1 (6,'Input pz' ,100,-1.,1.,0.) CALL HBOOK1 (7,'Input theta', 90, 0., 180., 0.0) CALL HBOOK1 (8,'Input phi', 90, -180., 180., 0.0) call HBOOK2 (9,'theta v phi',90,0.,180.,90,-180.,180.,0.) CALL HBOOK1 (10,'Exit Distance', 100, 0., 10., 0.0) CALL HBOOK1 (14,'exit px',100,-1.,1.,0.) CALL HBOOK1 (15,'exit py ',100,-1.,1.,0.) CALL HBOOK1 (16,'exit pz' ,100,-1.,1.,0.) CALL HBOOK1 (17,'exit theta', 90, 0., 180., 0.0) CALL HBOOK1 (18,'exit phi', 90, -180., 180., 0.0) call HBOOK2 (19,'exit theta v phi',90,0.,180.,90,-180.,180.,0.) do C=1,NEVENT IEVENT=IEVENT+1 * print *,'X0(1)=',xp(C),' X0(2)=',yp(C) X0(1)=xps(C) ! input coordinate of the user defined ray X0(2)=yps(C) pi180=0.017453292 THETA(C)=THETA(C)*pi180 * CALL ranlux(rano,1) PHI(C)=PHI(C)*pi180 V0(1) = SIN(THETA(C))*COS(PHI(C)) V0(2) = SIN(THETA(C))*SIN(PHI(C)) V0(3) = COS(THETA(C))*-1. * print *,theta(C),phi(C) px=V0(1) py=V0(2) pz=V0(3) CALL HFILL (3,xps(c),yps(c),1.) * print *,'px=',px(C),' py=',py(C),' pz=',(pz(C)*-1) * V0(1)=px(C) ! input direction of the user defined ray * V0(2)=py(C) ! not necessarily a unit vector * V0(3)=(pz(C)*-1) ! (it is normalized inside WICO) CALL HFILL (4,px,0.,1.) CALL HFILL (5,py,0.,1.) CALL HFILL (6,pz,0.,1.) * if (acos(pz(C)).gt.1.570796) then * th=180. - asin((sin(acos(pz(C))))*1.0/1.6)/3.14159*180 * else * th=asin((sin(acos(pz(C))))*1.0/1.6)/3.14159*180 * endif CALL HFILL (7,THETA(C)/pi180,0.,1.) CALL HFILL (8,PHI(C)/pi180,0.,1.) CALL HFILL (9,THETA(C)/pi180,PHI(C)/pi180,1.) C X0(1)=0. C X0(2)=0. * RMAX = TIN/2 * DO WHILE ((X0(1)**2+X0(2)**2) .GT. RMAX) * CALL ranlux(rano,1) * X0(1)=DBLE(RANO)*TIN-(TIN/2) ! input coordinate of the user defined ray * CALL ranlux(rano,1) * X0(2)=DBLE(RANO)*TIN-(TIN/2) * END DO * V0(1)=20. * V0(2)=20. C VMAX = 1.0875 C DO WHILE ((V0(1)**2+V0(2)**2) .GT. VMAX) * CALL ranlux(rano,1) * THETA= * COST = COS(rano*pi180*THETA)*(-1.0) * SINT2 = (1.-COST)*(1.+COST) * SINT = SQRT(SINT2) * CALL ranlux(rano,1) * PHI = 2*3.141592654*rano * SINP = SIN(PHI) * COSP = COS(PHI) C V0(1)=0. C V0(1)=DBLE(RANO)*2 - 1 ! input direction of the user defined ray C V0(2)=0. C V0(2)=DBLE(RANO)*2 - 1 C CAll ranlux(rano,1) C V0(3)=-1. C V0(3)=-DBLE(RANO) ! (it is normalized inside WICO) C END DO CALL WICO(X0,V0,X1,V1,D,N,IERR,NEVENT,IEVENT) if (IERR.eq.0) then nsuccess = nsuccess + 1 th1=180. -acos(V1(3))/pi180 c ph1=acos(V1(2)/sin(th1*pi180))/pi180 ph1=atan2(V1(1),V1(2))/pi180 px1=V1(1) py1=V1(2) pz1=V1(3) CALL HFILL (14,px1,0.,1.) CALL HFILL (15,py1,0.,1.) CALL HFILL (16,pz1,0.,1.) * if (acos(pz(C)).gt.1.570796) then * th=180. - asin((sin(acos(pz(C))))*1.0)/3.14159*180 * else * th=asin((sin(acos(pz(C))))*1.0)/3.14159*180 * endif CALL HFILL (17,TH1,0.,1.) CALL HFILL (18,ph1,0.,1.) CALL HFILL (19,th1,ph1,1.) endif CALL HISTFILL(X0,X1,IERR,D) IF(IERR.EQ.1) then PRINT*,'ray impact point out of entrance face: IERR=1' print *,th,X0(1),X0(2) endif * IF(IERR.EQ.2)PRINT*,'number of bounces greater than NMAX: IERR=2' * IF(IERR.EQ.3)PRINT*,' the ray comes back to the entrance face:IERR=3' * IF(IERR.EQ.4) then * PRINT*,' ray incident out of limit angle:IERR=4' * print *,V0(1),V0(2),V0(3),X0(1),X0(2),theta(C),phi(c) * endif * IF(IERR.EQ.0)THEN * PRINT*,' photon transmitted - IERR=0' * PRINT*,' output coordinate',X1(1),X1(2) * PRINT*,' output direction',V1(1),V1(2),V1(3) * PRINT*,' total path length',D * PRINT*,' total number of bounces',N * ENDIF C C=C+1 ENDDO c print *,' Number of successful Trials =',nsuccess print *,' Acceptance =', 1 float(nsuccess)/float(nevent) CALL HBOOKCLOSE STOP END SUBROUTINE HBOOKINIT C Create space in memory for histograms in HBOOK. C IMPLICIT NONE C Size of memory block reserved for histogram storage INTEGER HSTSIZ PARAMETER (HSTSIZ=500000) C C Define PAWC common INTEGER HSTBUF COMMON /PAWC/ HSTBUF(HSTSIZ) C CALL HLIMIT (HSTSIZ) RETURN END SUBROUTINE HBOOKSTART C Start up histograms C IMPLICIT NONE CALL HBOOK2 (1,'Input',100,-4.,4.,100,-4.,4.,0.) CALL HBOOK2 (2,'Output',100,-4.,4.,100,-4.,4.,0.) * CALL HBOOK2 (3,'Input x vs y',100,-3.,3.,100,-3.,3.,20.) RETURN END SUBROUTINE HISTFILL(X0,X1,IERR,D) C Fills the appropriate histograms with the appropriate values. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X0(2),X1(2) C PRINT*,' CHECK input coordinate',X0(1),X0(2) C PRINT*,' CHECK output coordinate',X1(1),X1(2) CALL HFILL (1,SNGL(X0(1)),SNGL(X0(2)),1.) IF(IERR.EQ.0)THEN CALL HFILL (2,SNGL(X1(1)),SNGL(X1(2)),1.) ENDIF CALL HFILL (10,SNGL(D),0.,1.) RETURN END SUBROUTINE HBOOKCLOSE C Saves the histograms created by HBOOK. C CALL HCDIR ('//PAWC',' ') CALL HRPUT (0,'wico.hbook','NT') RETURN END