program eta ***************************************************** * * this is a monte carlo program to simulate r+A--->eta+A, then * eta--->gamma1+e^- + e^+ , Sept 25, 2012, Liping * * ****************************************************************** double precision m_pi,m_pi0,m_p,m_eta,PI,m_he4,m_target,m_C12,m_e double precision P_r,P_p_a,P_p(3),theta_p,phi_p, P_g1_a, > P_g1(3),theta_g1,phi_g1 double precision P_pi0_a, P_pi0(3),theta_pi0,phi_pi0, > P_g2_a,P_g2(3),theta_g2,P_g3(3), > phi_g2, P_eta_a,P_eta(3),theta_eta,phi_eta double precision Pin_l,Pin_h,m(3),P(3),theta(3),phi(3),Px(3), > Py(3),Pz(3),Px_c(3),Py_c(3),Pz_c(3),philab(3), > theta_l(3),P_xy(3),ymax,sqm12,sqm23 double precision cos3_cc,E_pi0 double precision the_g3_c,P_g3_a, phi_g3_c,Px_g3_c,Py_g3_c, > Pz_g3_c,Px_g3,Py_g3,Pz_g3,theta_g3,phi_g3,P_xy_g3 double precision the_g4_c,P_g4_a, phi_g4_c,Px_g4_c,Py_g4_c, > Pz_g4_c,Px_g4,Py_g4,Pz_g4,theta_g4,phi_g4,P_xy_g4 double precision in_mass, Pxx,Pyy,Pzz,in_massr double precision an_max,dsig_dth real ebeam,cr_laget double precision x_g1,y_g1,x_g2,y_g2,x_g3,y_g3, > x_g4,y_g4,l_siz,up_siz integer pass_g1,pass_g2,pass_g3,pass_g4,N_pass,pass_sum integer loop,i,k integer eventID, ID_det double precision sig_E1,sig_p_g1,sig_E2,sig_p_g2, > sig_E3,sig_p_g3,sig_E4,sig_p_g4 real accept integer g1_out,g1_in,g1_pwo,g2_out,g2_in,g2_pwo,g3_out,g3_in, > g3_pwo,g4_out,g4_in,g4_pwo double precision x_g1_m,y_g1_m,x_g2_m,y_g2_m,x_g3_m,y_g3_m, > x_g4_m,y_g4_m double precision the_g1_m,phi_g1_m,p_g1_m,px_g1_m,py_g1_m,pz_g1_m double precision the_g2_m,phi_g2_m,p_g2_m,px_g2_m,py_g2_m,pz_g2_m double precision the_g3_m,phi_g3_m,p_g3_m,px_g3_m,py_g3_m,pz_g3_m double precision E_g2_m, E_g3_m,E_g2, E_g3 real ytest c real cross(1989,10),ytest,cross_t(4:12) integer ID_run, ID_out c integer Id_row,Id_Pr, ID_run, ID_out,flag1 double precision t,E_eta, E_p, mismass,P_p_m,E_p_m,the_p_m c dist_z (cm) is z distance from target to hycal, c siz_PWO (cm) is the half length of the PWO insertion c siz_cal (cm) is the half length of the calorimeter c siz_hol (cm) is the half length of the middle hole double precision dist_z,siz_pwo,siz_cal,siz_hol,z_vert integer N_shape PARAMETER (DEGRAD=0.01745329) PARAMETER (RADDEG=57.295779) real*4 ran_true PARAMETER (PI=3.141592653d0) c ini the particle mass in MeV data m_pi,m_p,m_pi0,m_eta,m_he4,m_C12,m_e /0.13956995d0, > 0.938272029d0,0.1349766d0,0.547853d0,3.73084d0,11.17788d0, > 0.5109989d-3/ call rluxgo(4,606968,0,0) c read in total events number to be generated print*,'input event number' read(5,*)n_ev write(6,*)'the total number of event is',n_ev c read beam energy print*,'input low limit of beam energy in GeV' read(5,*)Pin_l write(6,*)'low limit of beam energy is',Pin_l print*,'input high limit of beam energy in GeV' read(5,*)Pin_h write(6,*)'high limit of beam energy is',Pin_h print*,'input max theta angle for eta in Deg' read(5,*)an_max write(6,*)'max theta angle for eta is',an_max print*,'what kind of target? 1=H, 2=He4, 3=C12' read(5,*)k if(k.eq.1)then m_target=m_p elseif(k.eq.2)then m_target=m_he4 else m_target=m_C12 endif write(6,*)'target mass is',m_target c start to generate the event do loop=1,n_ev eventID=loop call ranlux(ran_true,1) P_r=dble(ran_true)*(Pin_h-Pin_l)+Pin_l c Id_Pr=idnint(P_r)-2 c***************************** C calculate the r+p--->eta+p according to cross section by Laget, Sept 10, 2011 ebeam=sngl(P_r) 100 call ranlux(ran_true,1) c Note: here angle is in deg theta_eta=dble(ran_true)*an_max call ranlux(ran_true,1) ytest=dble(ran_true*0.3) call CRSEC6(ebeam,theta_eta,cr_laget,t) dsig_dth=dble(cr_laget)*dsin(theta_eta/180.d0*PI) if(ytest.gt.dsig_dth) go to 100 c Note: convert angle into rad theta_eta=theta_eta/180.d0*PI c E_p=-t/(2.d0*m_target)+m_target c P_p_a=dsqrt(E_p**2-m_target**2) c E_eta=P_r+m_target-E_p c P_eta_a=dsqrt(E_eta**2-m_eta**2) c theta_p=dacos((P_r-P_eta_a*dcos(theta_eta))/P_p_a) call kin2_v(P_r,P_r,m_eta,m_target,theta_eta,E_eta) P_eta_a=dsqrt(E_eta**2-m_eta**2) E_p=P_r+m_target-E_eta P_p_a=dsqrt(E_p**2-m_target**2) theta_p=dacos((P_r-P_eta_a*dcos(theta_eta))/P_p_a) c************************** call ranlux(ran_true,1) phi_eta=dble(ran_true)*2.d0*PI if(phi_eta.lt.PI)then phi_p=phi_eta+PI else phi_p=phi_eta-PI endif P_eta(1)=P_eta_a*dsin(theta_eta)*dcos(phi_eta) P_eta(2)=P_eta_a*dsin(theta_eta)*dsin(phi_eta) P_eta(3)=P_eta_a*dcos(theta_eta) P_p(1)=P_p_a*dsin(theta_p)*dcos(phi_p) P_p(2)=P_p_a*dsin(theta_p)*dsin(phi_p) P_p(3)=P_p_a*dcos(theta_p) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC c start eta--->gamma+ e^- + e^+ c Note: g1 is gamma; g2 is e^-; g3 is e^+ c c m(1)=m_pi c m(2)=m_pi c m(3)=m_pi0 m(1)=m_pi0 m(2)=m_pi0 m(3)=m_pi0 C ymax is a number bigger than the max of P3_c*P1_cc c ymax=0.032 ymax=0.03 c subroutine body3(ymax,P0,M0,m,P,theta,phi) call body3(ymax,P_eta_a,m_eta,m,P,theta,phi) do i=1,3 Px_c(i)=P(i)*dsin(theta(i))*dcos(phi(i)) Py_c(i)=P(i)*dsin(theta(i))*dsin(phi(i)) Pz_c(i)=P(i)*dcos(theta(i)) Px(i)=Pz_c(i)*dsin(theta_eta)*dcos(phi_eta)+ > Px_c(i)*dcos(theta_eta)*dcos(phi_eta)- > Py_c(i)*dsin(phi_eta) Py(i)=Pz_c(i)*dsin(theta_eta)*dsin(phi_eta)+ > Px_c(i)*dcos(theta_eta)*dsin(phi_eta)+ > Py_c(i)*dcos(phi_eta) Pz(i)=Pz_c(i)*dcos(theta_eta)-Px_c(i)*dsin(theta_eta) if(P(i).gt.0.d0)then theta_l(i)=dacos(Pz(i)/P(i)) P_xy(i)=dsqrt(Px(i)**2+Py(i)**2) if(P_xy(i).gt.0.0d0) then philab(i)=dacos(Px(i)/P_xy(i)) if(Py(i).lt.0.0d0)philab(i)=2.d0*PI-philab(i) else call ranlux(ran_true,1) philab(i)=2.d0*PI*dble(ran_true) endif else call ranlux(ran_true,1) theta_l(i)=PI*dble(ran_true) call ranlux(ran_true,1) philab(i)=2.d0*PI*dble(ran_true) endif enddo c c c Note: output for eta--->gamma+ e^- + e^+ in the Lab frame c particle 1 is gamma, 2 is e^-, 3 is e^+ c c c c end of eta--->gamma+ e^- + e^+ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 206 continue c************** write in the generator output file open(unit=11, file='etadecay.dat') Id_run=9000 write(11,1)Id_run, loop 1 format(I5,1x,I8,1x,' 4') c write(11,10)m_eta c 10 format('1',' 17',1x, F7.5) c write(11,11)P_eta(1), P_eta(2), P_eta(3), E_eta c 11 format(' 0',1x, F10.6,1x,F10.6,1x,F10.6, 1x, F10.6) write(11,2)m_p 2 format('1',' 14',1x, F7.5) write(11,3)P_p(1), P_p(2), P_p(3), E_p 3 format(' 1',1x, F10.6,1x,F10.6,1x,F10.6, 1x, F10.6) write(11,4)m_pi0 4 format('2', ' 7', 1x, F7.5) write(11,5)Px(1), Py(1), Pz(1), P(1) 5 format(' 0',1x, F10.6,1x,F10.6,1x,F10.6, 1x, F10.6) write(11,8)m_pi0 8 format('3', ' 7', 1x, F7.5) write(11,9)Px(2), Py(2), Pz(2), P(2) 9 format(' 0',1x, F10.6,1x,F10.6,1x,F10.6, 1x, F10.6) write(11,12)0.0 12 format('4', ' 1', 1x, F7.5) write(11,13)Px(3), Py(3), Pz(3), P(3) 13 format(' 0',1x, F10.6,1x,F10.6,1x,F10.6, 1x, F10.6) enddo accept=real(N_pass)/real(n_ev) print*,'acceptance is ',accept print*,'N_pass is ',N_pass,', n_ev is ',n_ev if(ID_out.eq.1) then close(unit=11) endif write(6,*)'run is completed' stop end