C*****This subroutine fills the common NeutronsGR: Number_N, Neut_En(3), C teta_n(3), fi_n(3), C*****where C*****Number_N - number of neutrons in case of Gigantic resonance C (Number_N=0, if no); C*****NeutE_n(i) - energy of each generated neutron, MeV C*****teta_n(i) - teta (polar) angle in reference coordinate system C*****fi_n(i) - fi (azimuth) angle in reference coordinate system Input C parameters from GEANT: Step - GEANT step, cm En-gamma - gamma C energy in tungsten, GeV px_gamma, py_gamma, pz_gamma - gamma C momentum projections, GeV subroutine GPGRES(Step, En_g, px_g, py_g, pz_g) integer i real En_exc real cos_teta, sin_teta, cos_fi, sin_fi, En_n, alfa_n, beta_n C Input parameters * Energy of gamma, in GeV real En_gamma * Gamma momentum projections, multiplied by c, in GeV real px_gamma, py_gamma, pz_gamma * GEANT step, in cm real Step common/bind/En_bind(3) #include "gpgneut.inc" c * * Let's trasnform energy and momentum in MeV. Remember, GEANT is GeV * * open(77,file='testgig.dat',access='append') En_gamma=En_g/1000. px_gamma=px_g/1000. py_gamma=py_g/1000. pz_gamma=pz_g/1000. * * Exciting energy of nucleus in Gigantic resonance, MeV En_exc=En_gamma Number_N=0 do i=1,3 Neut_En(i)=0. teta_n(i)=0. fi_n(i)=0. enddo C Simulate number of neutrons in Gygantic resonance process C (Number_n=0 if no interaction) call GRESsimulate(En_gamma, Step, Number_N) if (Number_N.gt.0) then cos_teta=pz_gamma/sqrt(px_gamma**2+py_gamma**2+pz_gamma**2) sin_teta=sqrt(1.-cos_teta**2) cos_fi=px_gamma/sqrt(px_gamma**2+py_gamma**2) sin_fi=py_gamma/sqrt(px_gamma**2+py_gamma**2) do i=1, Number_N if (En_exc.gt.En_bind(i)) then call SpectrumSpace(Number_N, i, En_exc, & cos_teta, sin_teta, cos_fi, sin_fi, & En_n, alfa_n, beta_n) Neut_En(i)=En_n teta_n(i)=alfa_n fi_n(i)=beta_n En_exc=En_exc-En_bind(i)-En_n endif enddo endif * write(77,*) Number_N, step, en_gamma * close(77) RETURN end c************************************************************************ ************************************************************************* C*****Simulate generated number of C*****neutrons C****************************** subroutine GRESsimulate(En_gamma, Step, Number_N) c Output parameter: Number_N (as default Number_N=0) real En_gamma, Step integer Number_N real Sigma real*4 tmp real Lambda1, Lambda2, Lambda3 #include "gpgsigma.inc" if ((En_gamma.gt.9.097).and.(En_gamma.lt.28.608)) then CCCCC Calculate Lambda1 Sigma=-1. do i=1,63 if ((En_gamma.ge.energy1(i)).and.(En_gamma.lt.energy1(i+1))) & Sigma=cross1(i)+(cross1(i+1)-cross1(i))* & (En_gamma-energy1(i))/ & (energy1(i+1)-energy1(i)) enddo if (Sigma.le.0.) then Lambda1=1000. else 1 call GRNDM(tmp,1) if ((tmp.lt.1.E-8).or.(tmp.gt.(1.-1.E-8))) goto 1 Lambda1=-log(1.-tmp)/(Sigma*1.E-27*6.31*1.E22) endif CCCCCCCalculate Lambda2 Sigma=-1. do i=1,61 if ((En_gamma.ge.energy2(i)).and.(En_gamma.lt.energy2(i+1))) & Sigma=cross2(i)+(cross2(i+1)-cross2(i))* & (En_gamma-energy2(i))/ & (energy2(i+1)-energy2(i)) enddo if (Sigma.le.0.) then Lambda2=1000. else 2 call GRNDM(tmp,1) if ((tmp.lt.1.E-8).or.(tmp.gt.(1.-1.E-8))) goto 2 Lambda2=-log(1.-tmp)/(Sigma*1.E-27*6.31*1.E22) endif CCCCCCCalculate Lambda3 Sigma=-1. do i=1,26 if ((En_gamma.ge.energy3(i)).and.(En_gamma.lt.energy3(i+1))) & Sigma=cross3(i)+(cross3(i+1)-cross3(i))* & (En_gamma-energy3(i))/ & (energy3(i+1)-energy3(i)) enddo if (Sigma.le.0.) then Lambda3=1000. else 3 call GRNDM(tmp,1) if ((tmp.lt.1.E-8).or.(tmp.gt.(1.-1.E-8))) goto 3 Lambda3=-log(1.-tmp)/(Sigma*1.E-27*6.31*1.E22) endif if (Lambda1.eq.MIN(Lambda1, Lambda2, Lambda3, Step)) Number_N=1 if (Lambda2.eq.MIN(Lambda1, Lambda2, Lambda3, Step)) Number_N=2 if (Lambda3.eq.MIN(Lambda1, Lambda2, Lambda3, Step)) Number_N=3 endif RETURN end c************************************************************************ c************************************************************************ ******Simulate neutron energy and flight direction*********************** subroutine SpectrumSpace(Number_N, Num, En_exc, & cos_teta, sin_teta, cos_fi, sin_fi, En_n, alfa_n, beta_n) c Output parameters: En_n, alfa_n, beta_n common/bind/En_bind(3) real tmp(2) real Pi, Temperature, lev_dens, cos_alfa, sin_alfa, beta real cos_alfa_n, cos_beta_n, sin_beta_n integer Number_N, Num real En_exc real cos_teta, sin_teta, cos_fi, sin_fi, En_n, alfa_n, beta_n Pi=4.*atan(1.) if (Number_N.eq.1) then En_n=En_exc-En_bind(1) else if ((Number_N.eq.2).and.(Num.eq.1)) then C************************************************************************ C SIMULATION OF NEUTRON EVAPORATION SPECTRUM The reference book C is: V.C. Barashenkov, V.D. Toneev, "The interactions of high C energy particles and atomic nuclei with nuclei", in Russian C Calculation of density levels parameter (see (6.31) from reference C book) lev_dens=0.085*184. C Calculation of nucleus temperature (see (6.30) from reference book) C (see (6.64) and comments after) Temperature=sqrt((En_exc-En_bind(Num))/lev_dens) C Neutron evaporation energy spectrum simulation by exception method 4 call GRNDM(tmp,2) C The average binding energy of neutron x_tmp=(En_exc-En_bind(1)-En_bind(2))*tmp(1) y_tmp=(exp(-1.)/Temperature)*tmp(2) if (y_tmp.gt.((x_tmp/Temperature**2)*exp(-x_tmp/Temperature))) & goto 4 En_n=x_tmp else if ((Number_N.eq.2).and.(Num.eq.2)) then En_n=En_exc-En_bind(2) else if ((Number_N.eq.3).and.(Num.eq.1)) then C************************************************************************ C SIMULATION OF NEUTRON EVAPORATION SPECTRUM The reference book is: C V.C. Barashenkov, V.D. Toneev, "The interactions of high energy C particles and atomic nuclei with nuclei", in Russian Calculation C of density levels parameter (see (6.31) from reference book) lev_dens=0.085*184. C Calculation of nucleus temperature (see (6.30) from reference book) C (see (6.64) and comments after) Temperature=sqrt((En_exc-En_bind(Num))/lev_dens) C Neutron evaporation energy spectrum simulation by exception method 5 call GRNDM(tmp,2) C Average binding energy of neutron x_tmp=(En_exc-En_bind(1)-En_bind(2)-En_bind(3))*tmp(1) y_tmp=(exp(-1.)/Temperature)*tmp(2) if (y_tmp.gt.((x_tmp/Temperature**2)*exp(-x_tmp/Temperature))) & goto 5 En_n=x_tmp else if ((Number_N.eq.3).and.(Num.eq.2)) then C************************************************************************ C SIMULATION OF NEUTRON EVAPORATION SPECTRUM The reference book C is: V.C. Barashenkov, V.D. Toneev, "The interactions of high C energy particles and atomic nuclei with nuclei", in Russian C Calculation of density levels parameter (see (6.31) from reference C book) lev_dens=0.085*183. C Calculation of nucleus temperature (see (6.30) from reference book) C (see (6.64) and comments after) Temperature=sqrt((En_exc-En_bind(Num))/lev_dens) C Neutron evaporation energy spectrum simulation by exception method 6 call GRNDM(tmp,2) C Average binding energy of neutron x_tmp=(En_exc-En_bind(2)-En_bind(3))*tmp(1) y_tmp=(exp(-1.)/Temperature)*tmp(2) if (y_tmp.gt.((x_tmp/Temperature**2)*exp(-x_tmp/Temperature))) & goto 6 En_n=x_tmp else if ((Number_N.eq.3).and.(Num.eq.3)) then En_n=En_exc-En_bind(3) endif C************************************************************************ C Neutron direction simulation in gamma coordinate system call GRNDM(tmp,2) cos_alfa=-1.+2.*tmp(1) beta=2.*Pi*tmp(2) sin_alfa=sqrt(1.-cos_alfa**2) C Return to the GEANT reference coordinate system cos_alfa_n=-sin_alfa*cos(beta)*(sin_teta)+cos_alfa*(cos_teta) cos_beta_n=sin_alfa*cos(beta)*(cos_teta*cos_fi)- & sin_alfa*sin(beta)*(sin_fi)+cos_alfa*(sin_teta*cos_fi) sin_beta_n=sin_alfa*cos(beta)*(cos_teta*sin_fi)+ & sin_alfa*sin(beta)*(cos_fi)+cos_alfa*(sin_teta*sin_fi) if (abs(cos_beta_n).lt.1.E-6) cos_beta_n=0. if (abs(sin_beta_n).lt.1.E-6) sin_beta_n=0. if ((cos_beta_n.eq.0.).and.(sin_beta_n.eq.0.)) beta_n=0. alfa_n=acos(cos_alfa_n) if ((cos_beta_n.ge.0.).and.(sin_beta_n.gt.0.)) & beta_n=atan(sin_beta_n/cos_beta_n) if ((cos_beta_n.ge.0.).and.(sin_beta_n.lt.0.)) & beta_n=2.*Pi+atan(sin_beta_n/cos_beta_n) if ((cos_beta_n.lt.0.).and.(sin_beta_n.le.0.)) & beta_n=Pi+atan(sin_beta_n/cos_beta_n) if ((cos_beta_n.lt.0.).and.(sin_beta_n.ge.0.)) & beta_n=Pi+atan(sin_beta_n/cos_beta_n) RETURN END