SUBROUTINE GUKINE C Generate the kinematics -- store tracks and vertices for C this event. IMPLICIT NONE # include "uglpar.inc" # include "geant321/gclist.inc" # include "geant321/gcunit.inc" # include "geant321/gcnum.inc" # include "ucglob.inc" # include "uckine.inc" # include "ucgeio.inc" # include "ucvert.inc" # include "ucconf.inc" # include "ucphys.inc" # include "ucntup.inc" # include "geant321/gcvolu.inc" # include "geant321/gccuts.inc" # include "geant321/gcking.inc" # include "geant321/gctrak.inc" # include "geant321/gcflag.inc" # include "geant321/gcnum.inc" c added by Jeff for access to mode # include "bfield.inc" INTEGER NVERT,IP,NTT,I_PARTICLE LOGICAL GOOD_EVENT INTEGER NUBUF,IADR REAL UBUF(3) REAL XNORM,YNORM,RANDOM,VTX(3),PLB(3),P_LAB(5),GEOMWT REAL BETA(3),SIGMA_DIPOLE REAL SIGMA_PI,THETA_PI,SIGMA_EEP,THETA_N,PHI,THETA REAL PT2,PT2_X,PT2_Y,PT2_Z,T_PROTON1,T_PROTON2 REAL PN_LAB,PN_Y,PN_X,PN_Z REAL EP1_CM,EP2_CM,Q_VAL,TOT_KIN REAL E_THETA,E_PHI,E_MOMENTUM REAL THETA_DEG,PHI_DEG,D_LIMIT,U_LIMIT REAL TH_CM,E_ENERGY,THETA_M1,MOL_CROSS REAL E_PRIME REAL PHI_M1 DOUBLE PRECISION THETA_N_D,W_N_D,E_L_D REAL PR1_PX_LAB,PR1_PY_LAB, PR1_PZ_LAB,PR1_P_LAB,PR1_E_LAB REAL PR2_PX_LAB,PR2_PY_LAB, PR2_PZ_LAB,PR2_P_LAB,PR2_E_LAB REAL AN_PX_LAB,AN_PY_LAB, AN_PZ_LAB,AN_P_LAB,AN_E_LAB REAL THETA_PR1,PHI_PR1,THETA_PR2,PHI_PR2 REAL T_TWO_PROTONS,Q_REACTION,THETA_PRN,PHI_PRN REAL C_EN,V_EN_X,V_EN_Y,V_EN_Z,V_EN,THETA_V_EN,PHI_V_EN REAL V_X,V_Y,V_Z,PHI_EL,E_X_LAB,E_Y_LAB,E_Z_LAB REAL CS_EN,OO2,ELECTRON_E_LAB,ELECTRON_P_LAB,THETA_EL REAL ELECTRON_PX_LAB,ELECTRON_PY_LAB,ELECTRON_PZ_LAB C Variables for GFPART (Extracts the constants describing a particle) REAL CHARGE,TLIFE,UB,ENERGY INTEGER ITRTYPE,NWB CHARACTER*20 NAIPART C Mass used in this subroutine REAL MPROTON,MNEUTRON,MELEC,MPIP,MPIM,MPI0,MDEUTERON,MHE3 c momentum range for inclusive reaction real p_min,p_max c extra temporary variables for Brems generator real p_nu,theta_cm,ppi,thetapi integer iso c variables used in calculating bremsstrahlung corrections real incident_electron_energy_corrected logical lstrak_(6) real unif real chi real alphaQED,PI_ parameter (alphaQED=7.2974e-3,PI_=3.1415926) ** ** For Short Range Correlations: ** REAL FSPACE(200), FERMI_DIST,FINT REAL PT1_X,PT1_Y,PT1_Z,PT1 REAL XLOW,XHIGH,T_NEUTRON,T_PROTON,RMOM REAL GAMMA,EP_CM,EN_CM,PHI_N EXTERNAL FERMI_DIST COMMON / FUNINT/ FINT REAL XXMOTT EXTERNAL XXMOTT REAL CROSS_MOTT,G_D_v,TAU,FAC,T2THE REAL G_P_E,S_n,G_P_M,G_N_M,SIGMA_DIPOLE_N,G_N_EL C Initialize masses CALL GFPART(3,NAIPART,ITRTYPE,MELEC,CHARGE,TLIFE,UB,NWB) CALL GFPART(13,NAIPART,ITRTYPE,MNEUTRON,CHARGE,TLIFE,UB,NWB) CALL GFPART(14,NAIPART,ITRTYPE,MPROTON,CHARGE,TLIFE,UB,NWB) CALL GFPART(7,NAIPART,ITRTYPE,MPI0,CHARGE,TLIFE,UB,NWB) CALL GFPART(8,NAIPART,ITRTYPE,MPIP,CHARGE,TLIFE,UB,NWB) CALL GFPART(9,NAIPART,ITRTYPE,MPIM,CHARGE,TLIFE,UB,NWB) CALL GFPART(45,NAIPART,ITRTYPE,MDEUTERON,CHARGE,TLIFE,UB,NWB) CALL GFPART(49,NAIPART,ITRTYPE,MHE3,CHARGE,TLIFE,UB,NWB) C Initialize the HITS flag of the detector CALL VZERO(HDETEC,MAXDET*MAXEL) C Initialize the list of detector's order... CALL VZERO(P_ORDER,MAXPAR) C Calculate and store vertex (X and Y) IF (BSIGX.GE.0..AND.BSIGY.GE.0.) THEN CALL GRANOR (XNORM,YNORM) VTX(1) = BSIGX*XNORM+BPOSX VTX(2) = BSIGY*YNORM+BPOSY ELSE VTX(1) = -BSIGX*2.*(RANDOM()-0.5)+BPOSX VTX(2) = -BSIGY*2.*(RANDOM()-0.5)+BPOSY ENDIF TAR_X = VTX(1) TAR_Y = VTX(2) RE_X = VTX(1) RE_Y = VTX(2) C Position vertex along Z according to user selection IF (IPTY.EQ.0) THEN C Randomly position the vertex Z within the target. VTX(3) = 2.*(RANDOM()-0.5)*BXTARG(3) ELSE C Position the vertex at the one of the ends or in C the middle of the target. VTX(3) = MIN(1,MAX(-1,IPTY-2))*BXTARG(3) ENDIF VTX(3) = VTX(3)+Z_TAR TAR_Z = VTX(3) RE_Z = VTX(3) C To determine the energy of the electron incident on the scattering C vertex inside the target, the energy loss up to the scattering C event (inside the target material and walls) must be simulated. C The way this is done here is to first determine the scattering vertex C position, then track the beam electron in reverse towards the source. C Any secondaries produced are going the wrong way, so secondaries are C disabled during this pre-event tracking step. C August 30, 2004: richard.t.jones@uconn.edu call GSVERT (VTX,0,0,0,0,NVERT) I_PARTICLE = 3 P_LAB(1) = 0 P_LAB(2) = 0 P_LAB(3) = -EBEAM call GSKINE (P_LAB,I_PARTICLE,NVERT,0,0,NTT) do IP=1,6 lstrak_(IP) = LSTRAK(IP) LSTRAK(IP) = .false. enddo call GTREVE do IP=1,6 LSTRAK(IP) = lstrak_(IP) enddo incident_electron_energy_corrected = GETOT NTRACK = 0 NVERTX = 0 IEOTRI = 0 TOFG = 0 C **************************** C * Event from an event file * C **************************** IF (I_REACTION.EQ.0) THEN C Calculate and store vertex CALL GSVERT (VTX,0,0,0,0,NVERT) IF (NVERT.EQ.0) GO TO 230 READ(IUKINE,*)NP,Q__2 PRINT*,'nombre de particule ',NP DO IP=1,NP READ(IUKINE,*)I_PARTICLE,P_LAB(1),P_LAB(2),P_LAB(3),weight_n CALL GFPART(I_PARTICLE,NAIPART,ITRTYPE,AMASS(1),CHARGE, $ TLIFE,UB,NWB) P_L_X(I_PARTICLE)=P_LAB(1) P_L_Y(I_PARTICLE)=P_LAB(2) P_L_Z(I_PARTICLE)=P_LAB(3) P_L_TOT(I_PARTICLE)=SQRT(P_LAB(1)**2+P_LAB(2)**2+ $ P_LAB(3)**2) P_LAB(5)=P_L_TOT(I_PARTICLE) E_L(I_PARTICLE)=SQRT(P_LAB(5)**2+AMASS(1)*2) c WEIGHT_N=1 IF(LPTRAK(I_PARTICLE))THEN PLB(1)=P_L_X(I_PARTICLE) PLB(2)=P_L_Y(I_PARTICLE) PLB(3)=P_L_Z(I_PARTICLE) CALL GSKINE (PLB,I_PARTICLE,NVERT,0,0,NTT) PRINT*,'particule ',I_PARTICLE,' prise ? ',NTT IF (NTT.EQ.0) GO TO 220 NUBUF=2 UBUF(1)=0 UBUF(2)=1 IADR=0 CALL GSKINU(NTT,NUBUF,UBUF,IADR) ELSE WRITE (CHMAIL,1) I_PARTICLE 1 FORMAT (' WARNING : Particle ',I3, $ ' defined in the event file and not tracked !') CALL GMAIL(0,0) ENDIF ENDDO IF (NTRACK.GT.0) THEN I_EVENT=I_EVENT+1 PRINT*,'EVENT ',I_EVENT ENDIF RETURN ENDIF IF (I_REACTION.GT.10) GOTO 320 GO TO (10,20,70,120,170,601,701,801) I_REACTION C **************************************** C * This part for electromagnetic shower * C **************************************** C PROBLEM: MANY items in ntuple are not defined for this C reaction! 10 continue C Calculate and store vertex : C Electron's origin is just before the target C VTX(3) = -75.0 VTX(1)=CER_X(I_EVENT) VTX(2)=CER_Y(I_EVENT) CALL GSVERT (VTX,0,0,0,0,NVERT) IF (NVERT.EQ.0) GO TO 230 C C Calculate and store kinematics C Electron: I_PARTICLE = 3 C Define electron kinematics to be that of the incident beam. CCERENKOV PHI=RANDOM()*ACOS(-1.)*2 PHI_DEG=CER_PHI(I_EVENT) c IF(PHI_DEG.LT.0.) PHI_DEG=PHI_DEG+180. PHI_DEG=PHI_DEG+90. PHI=PHI_DEG*3.1415927/180. WEIGHT_N=CER_W(I_EVENT) IF(BDIV.LT.0) THEN THETA=SQRT(-LOG(RANDOM())/LOG(2.))*ABS(BDIV)*ACOS(-1.)/180. ELSE THETA_DEG=CER_THE(I_EVENT) THETA=THETA_DEG*3.1415927/180. ENDIF c print *, VTX(1),VTX(2),VTX(3),PHI_DEG,THETA_DEG IF (EBEAM.EQ.3.AND.IPTY.EQ.3) THEN C Take into account energy loss inside the target (using C "Maurice's NIKHEF program) : 2989.9 += 25.59 MeV => C 3 GeV - 2.965 GeV ENERGY = 3-RANDOM()*0.035 ELSE ENERGY = EBEAM ENDIF P_LAB(1) = ENERGY*SIN(THETA)*COS(PHI) P_LAB(2) = ENERGY*SIN(THETA)*SIN(PHI) P_LAB(3) = ENERGY*COS(THETA) P_LAB(4) = ENERGY P_LAB(5) = ENERGY C Select the direction of the beam through the geometry IF (IBDIR.EQ.0) THEN PLB(1) = P_LAB(1) PLB(3) = P_LAB(3) ELSE PLB(1) = -P_LAB(1) PLB(3) = -P_LAB(3) ENDIF PLB(2) = P_LAB(2) C Save momenta and enerties for possible ntuple output. P_L_X(I_PARTICLE) = PLB(1) P_L_Y(I_PARTICLE) = PLB(2) P_L_Z(I_PARTICLE) = PLB(3) E_L(I_PARTICLE) = P_LAB(4) P_L_TOT(I_PARTICLE) = P_LAB(5) I_EVENT = I_EVENT+1 IF (LPTRAK(I_PARTICLE)) THEN CALL GSKINE (PLB,I_PARTICLE,NVERT,0,0,NTT) IF (NTT.EQ.0) GO TO 220 NUBUF=2 UBUF(1)=0 UBUF(2)=1 IADR=0 CALL GSKINU(NTT,NUBUF,UBUF,IADR) ENDIF RETURN C C ****************************************************** C * This part for electron-proton and electron-neutron * C * elastic scattering * C ****************************************************** C 20 continue C Calculate and store vertex CALL GSVERT (VTX,0,0,0,0,NVERT) IF (NVERT.EQ.0) GO TO 230 C PRINT*,'NVERT NTRACK',NVERT,NTRACK C C Call GENBOD to get one event C Input C Number of outgoing particles NP = 2 C Electron mass AMASS(1) = MELEC C Electron energy at vertex E_ENERGY = incident_electron_energy_corrected C Proton mass AMASS(2) = MPROTON C Weight KGENEV = 1 C Total cm en. TECM = SQRT(AMASS(2)**2+2.*E_ENERGY*AMASS(2)) BETA(1) = 0. BETA(2) = 0. BETA(3) = -E_ENERGY/(E_ENERGY+AMASS(2)) C If NTT is set to 0., and the first event does not fit the cut, C GSKINE is not called and a "false error" is detected... C NTT = 0. NTT=1. I_EVENT = I_EVENT+1 IF (ISOTROPE.GT.0) THEN C This part is to have an isotropic distribution in the Lab. C Which is necessary to compute rates... C Each theta is between 0 and PI ; each phi is between 0 C and 2.PI C C Electron IF (ISOTROPE.EQ.1) THEN E_THETA=RANDOM()*(E_ANMAX-E_ANMIN)+E_ANMIN ELSE C C This section was added to enable a more realistic generation of C backgrounds for elastic scattering by putting the 1/theta_lab**3 C approximate weight factor from the elastic cross section into the C generator. The departure of cross section from 1/theta_lab**3 is C reflected in the weight. Bounds on the allowed values of theta_lab C and phi_lab are the same as for the ISOTROPE=1 case. C August 29, 2004: richard.t.jones@uconn.edu C unif=RANDOM() E_THETA=1/sqrt(unif/E_ANMAX**2+(1-unif)/E_ANMIN**2) ENDIF E_PHI=RANDOM()*(PHIMAX-PHIMIN)+PHIMIN E_MOMENTUM=E_ENERGY/(1+E_ENERGY/MPROTON*(1-COS(E_THETA))) Q__2=4*E_ENERGY*E_MOMENTUM*SIN(E_THETA/2)**2 C C Here we must account for internal bremsstrahlung at the elastic C scattering vertex. The correction to second-order in alphaQED is C C d sigma /d sigma\ / 2 alphaQED \ C ------- = |---------| | 1 - ---------- ln(E/Kcut) Chi(q2) | C d Omega \d Omega/ 0 \ pi / C C (see Bjorken and Drell "Relativistic Quantum Mechanics, Eq. 8.81) C where Kcut is the elastic infrared cutoff for bremsstrahlung and C Chi(q2) -> ln(-q2/m**2)-1 [ibid. Eq. 8.66] in the relativistic limit. C The following code generates energy loss consistent with this cross C section. Strictly speaking it is only valid for Kcut << E, but this C restriction holds for the situation of interest for Qweak simulations. C August 30, 2004: richard.t.jones@uconn.edu C unif=RANDOM() chi=log(Q__2/MELEC**2)-1 chi=max(chi, 0.00001) E_ENERGY=E_ENERGY*(1-exp(-unif*PI_/(2*alphaQED*chi))) C C To obtain consistent kinematics, we treat the internal bremsstrahlung C photon as if it were colinear with the outgoing electron, and just C attenuate the energy of the scattered electron. This leaves the C q2 of the vertex unchanged, and it corresponds to the peaking C approximation in the case of forward scattering, but away from the C forward direction it ignores small momentum components transverse C to the outgoing electron direction. At some point this should be C checked because it affects the experimental q2 acceptance. C August 30, 2004: richard.t.jones@uconn.edu C E_MOMENTUM=E_ENERGY/(1+E_ENERGY/MPROTON*(1-COS(E_THETA))) PLB(1)=E_MOMENTUM*COS(E_PHI)*SIN(E_THETA) PLB(2)=E_MOMENTUM*SIN(E_PHI)*SIN(E_THETA) PLB(3)=E_MOMENTUM*COS(E_THETA) C IF (LPTRAK(3)) THEN CALL GSKINE(PLB,3,NVERT,0,0,NTT) NUBUF=2 UBUF(1)=0 UBUF(2)=1 IADR=0 CALL GSKINU(NTT,NUBUF,UBUF,IADR) ENDIF C C Proton PLB(1)=-PLB(1) PLB(2)=-PLB(2) PLB(3)=E_ENERGY-PLB(3) C IF (LPTRAK(14)) THEN CALL GSKINE(PLB,14,NVERT,0,0,NTT) NUBUF=2 UBUF(1)=0 UBUF(2)=1 IADR=0 CALL GSKINU(NTT,NUBUF,UBUF,IADR) ENDIF C if(itarget.eq.1)then CALL ELASTIC_CROSS_SECTION_PROTON(E_THETA,WEIGHT_N) else call quasi_elastic(e_theta,weight_n,1,1) endif if (ISOTROPE.eq.1) then WEIGHT_N=WEIGHT_N*SIN(E_THETA) else WEIGHT_N=WEIGHT_N*SIN(E_THETA)*E_THETA**3* &(1/E_ANMIN**2-1/E_ANMAX**2)/(E_ANMAX-E_ANMIN)/2 endif ELSE 30 CALL GENBOD GOOD_EVENT=.FALSE. C Proton I_PARTICLE = 14 CALL UGPKIN(I_PARTICLE,2,BETA,P_LAB,PLB,THETA_P,GEOMWT) QQ__2 = 2.*0.93827**2*(SQRT(1.+(P_LAB(5)/0.93827)**2)-1.) IF (P_LAB(1).NE.0) THEN PHI_P=ATAN(P_LAB(2)/P_LAB(1)) IF (P_LAB(1).LT.0) PHI_P=PHI_P+3.14159 ELSE IF (P_LAB(2).GE.0) PHI_P=3.14159/2. IF (P_LAB(2).LT.0) PHI_P=-3.14159/2. ENDIF IF (LPTRAK(I_PARTICLE)) THEN IF ((THETA_P.GT.THETAMIN).AND.(THETA_P.LT.THETAMAX).AND. $ (QQ__2.GT.QMIN).AND.(QQ__2.LT.QMAX).AND. $ (PHI_P.GT.PHIMIN).AND.(PHI_P.LT.PHIMAX)) THEN CALL GSKINE(PLB,I_PARTICLE,NVERT,0,0,NTT) IF (NTT.EQ.0) GO TO 220 GOOD_EVENT=.TRUE. NUBUF=2 UBUF(1)=0 UBUF(2)=1 IADR=0 CALL GSKINU(NTT,NUBUF,UBUF,IADR) ENDIF C The above line was changed to (to quiet the compiler): IF (.NOT. GOOD_EVENT) GOTO 30 ENDIF C C Electron: I_PARTICLE = 3 THETA_E=E_THETA C CALL UGPKIN(I_PARTICLE,1,BETA,P_LAB,PLB,THETA_E,GEOMWT) C Find the weight if(itarget.eq.1)then CALL ELASTIC_CROSS_SECTION_PROTON(THETA_E,SIGMA_DIPOLE) else call quasi_elastic(theta_e,sigma_dipole,1,1) endif WEIGHT_N = SIGMA_DIPOLE/WT IF (LPTRAK(I_PARTICLE)) THEN GOOD_EVENT=.FALSE. IF ((THETA_E.GT.E_ANMIN).AND.(THETA_E.LT.E_ANMAX)) THEN CALL GSKINE(PLB,I_PARTICLE,NVERT,0,0,NTT) IF (NTT.EQ.0) GO TO 220 NUBUF=2 UBUF(1)=0 UBUF(2)=1 IADR=0 CALL GSKINU(NTT,NUBUF,UBUF,IADR) GOOD_EVENT=.TRUE. ENDIF ENDIF C Normally, return after elastic event generation. But C in special plot mode, go on to generate an inelastic event C for the plot. If EITHER electron or proton did not satify C the cut, generate a new one IF(.NOT.GOOD_EVENT)THEN C Remove any other particule whose cut was satisfied NTRACK=0 GO TO 30 ELSE Q__2=QQ__2 ENDIF ENDIF RETURN C C ****************************************************** C * This part for electron-proton and electron-neutron * C * inelastic scattering * C * e+p ==> e+n+(pi+) * C * Modified by J. Mammei 1-28-2008 * C ****************************************************** C 70 continue C Calculate and store vertex CALL GSVERT (VTX,0,0,0,0,NVERT) IF (NVERT.EQ.0) GO TO 230 C PRINT*,'NVERT NTRACK',NVERT,NTRACK C C Call GENBOD to get one event C Input C Number of outgoing particles NP = 2 C Incident Electron mass AMASS(1) = MELEC C Electron energy at vertex E_ENERGY = incident_electron_energy_corrected C Proton mass AMASS(2) = MPROTON C Pion mass AMASS(3) = MPIP C Weight KGENEV = 1 C Total cm en. TECM = SQRT(AMASS(1)**2+AMASS(2)**2+2.*E_ENERGY*AMASS(2)) BETA(1) = 0. BETA(2) = 0. BETA(3) = -SQRT(E_ENERGY**2-AMASS(1)**2)/(E_ENERGY+AMASS(2)) C If NTT is set to 0., and the first event does not fit the cut, C GSKINE is not called and a "false error" is detected... C NTT = 0. NTT=1. I_EVENT = I_EVENT+1 C Electron I_PARTICLE = 3 IF (LPTRAK(I_PARTICLE)) THEN C Throw flat in theta and phi, as we do for elastics and mollers E_THETA=RANDOM()*(E_ANMAX-E_ANMIN)+E_ANMIN E_PHI=RANDOM()*(PHIMAX-PHIMIN)+PHIMIN C Also throw flat in E'; define E'_min and E'_max here temporarily C later should put into bash script? C EP_MIN = MELEC C EP_MAX = EBEAM C EP_MAX = E_ENERGY-AMASS(3)*(1.0+0.5*AMASS(3)/AMASS(2)) C EP_MAX = EP_MAX/(1.0+(2.0*E_ENERGY/AMASS(2))*SIN(E_THETA/2)**2) E_PRIME=RANDOM()*(EBEAM-MELEC)+MELEC E_MOMENTUM=SQRT(E_PRIME**2-AMASS(1)**2) PLB(1)=E_MOMENTUM*COS(E_PHI)*SIN(E_THETA) PLB(2)=E_MOMENTUM*SIN(E_PHI)*SIN(E_THETA) PLB(3)=E_MOMENTUM*COS(E_THETA) C Do we need this? C CALL UGPKIN(I_PARTICLE,1,BETA,P_LAB,PLB,THETA_E,GEOMWT) C Inelastic cross section CALL SIGMA_EEPRIME (E_ENERGY,E_PRIME,E_THETA,SIGMA_EEP) EL_WEIGHT = SIGMA_EEP WEIGHT_N=EL_WEIGHT*SIN(E_THETA)*(EBEAM-MELEC) Q__2=4*E_ENERGY*E_PRIME*SIN(E_THETA/2)**2 C IF (LPTRAK(I_PARTICLE)) THEN CALL GSKINE(PLB,3,NVERT,0,0,NTT) NUBUF=2 UBUF(1)=0 UBUF(2)=1 IADR=0 CALL GSKINU(NTT,NUBUF,UBUF,IADR) ENDIF ENDIF RETURN C ****************************************************** C * This part for electron-proton and electron-neutron * C * inelastic scattering * C * e+p ==> e+p+(pi+)+(pi-) * C ****************************************************** C 120 continue C Calculate and store vertex CALL GSVERT (VTX,0,0,0,0,NVERT) IF (NVERT.EQ.0) GO TO 230 C Call GENBOD to get one event C Input: C Number of outgoing particles NP = 4 C Electron mass AMASS(1) = MELEC C Proton mass AMASS(2) = MPROTON C Pion mass AMASS(3) = MPIP C Pion mass AMASS(4) = MPIM C Weight KGENEV = 1 C Total cm en. TECM = SQRT(AMASS(2)**2+2.*EBEAM*AMASS(2)) BETA(1) = 0. BETA(2) = 0. BETA(3) = -EBEAM/(EBEAM+AMASS(2)) I_EVENT = I_EVENT+1 130 CALL GENBOD C Check on Electron: I_PARTICLE = 3 CALL UGPKIN(I_PARTICLE,1,BETA,P_LAB,PLB,THETA_E,GEOMWT) IF ((THETA_E.LT.E_ANMIN).OR.(THETA_E.GT.E_ANMAX)) GOTO 130 E_PHI=ACOS(P_LAB(1)/SQRT(P_LAB(1)**2+P_LAB(2)**2)) IF(P_LAB(2).LT.0.) E_PHI=6.2831853-E_PHI IF ((E_PHI.LT.6.2831853+PHIMIN).AND.(E_PHI.GT.PHIMAX)) GOTO 130 C C Negative pion: I_PARTICLE = 9 IF (LPTRAK(I_PARTICLE)) THEN CALL UGPKIN(I_PARTICLE,3,BETA,P_LAB,PLB,THETA_PI,GEOMWT) IF ((THETA_PI.LT.E_ANMIN).OR.(THETA_PI.GT.E_ANMAX)) GOTO 130 C Find the weight THETA_PI = ACOS(P_LAB(3)/P_LAB(5)) IP = -2 CALL EPC_O (IP,EBEAM,P_LAB(5),THETA_PI,SIGMA_PI) C ns EL_WEIGHT = SIGMA_PI/TOT_MOTT_CROSS/GEOMWT EL_WEIGHT = SIGMA_PI WEIGHT_N=EL_WEIGHT CALL GSKINE(PLB,I_PARTICLE,NVERT,0,0,NTT) NUBUF=2 UBUF(1)=0 UBUF(2)=1 IADR=0 CALL GSKINU(NTT,NUBUF,UBUF,IADR) IF (NTT.EQ.0) GO TO 220 ENDIF C C Positive pion: I_PARTICLE = 8 IF (LPTRAK(I_PARTICLE)) THEN CALL UGPKIN(I_PARTICLE,4,BETA,P_LAB,PLB,THETA_PI,GEOMWT) CALL GSKINE(PLB,I_PARTICLE,NVERT,0,0,NTT) NUBUF=2 UBUF(1)=0 UBUF(2)=1 IADR=0 CALL GSKINU(NTT,NUBUF,UBUF,IADR) IF (NTT.EQ.0) GO TO 220 ENDIF C C Proton: I_PARTICLE = 14 IF (LPTRAK(I_PARTICLE)) THEN CALL UGPKIN(I_PARTICLE,2,BETA,P_LAB,PLB,THETA_P,GEOMWT) CALL GSKINE(PLB,I_PARTICLE,NVERT,0,0,NTT) NUBUF=2 UBUF(1)=0 UBUF(2)=1 IADR=0 CALL GSKINU(NTT,NUBUF,UBUF,IADR) IF (NTT.EQ.0) GO TO 220 ENDIF C C Electron: I_PARTICLE = 3 IF (LPTRAK(I_PARTICLE)) THEN CALL UGPKIN(I_PARTICLE,1,BETA,P_LAB,PLB,THETA_E,GEOMWT) CALL GSKINE(PLB,I_PARTICLE,NVERT,0,0,NTT) NUBUF=2 UBUF(1)=0 UBUF(2)=1 IADR=0 CALL GSKINU(NTT,NUBUF,UBUF,IADR) IF (NTT.EQ.0) GO TO 220 ENDIF RETURN C C ****************************************************** C * This part for electron-proton and electron-neutron * C * inelastic scattering * C * e+p ==> e+p+(pi0) * C ****************************************************** C 170 continue C Calculate and store vertex CALL GSVERT (VTX,0,0,0,0,NVERT) IF (NVERT.EQ.0) GO TO 230 C C Call GENBOD to get one event C Input: C Number of outgoing particles NP = 3 C Electron mass AMASS(1) = MELEC C Proton mass AMASS(2) = MPROTON C Pion mass AMASS(3) = MPI0 C Weight KGENEV = 1 C Total cm en. TECM = SQRT(AMASS(2)**2+2.*EBEAM*AMASS(2)) BETA(1) = 0. BETA(2) = 0. BETA(3) = -EBEAM/(EBEAM+AMASS(2)) I_EVENT = I_EVENT+1 180 CALL GENBOD C C Check on PROTON: I_PARTICLE = 14 CALL UGPKIN(I_PARTICLE,2,BETA,P_LAB,PLB,THETA_P,GEOMWT) IF ((THETA_P.LT.THETAMIN).OR.(THETA_P.GT.THETAMAX)) GOTO 180 PHI_P=ACOS(P_LAB(1)/SQRT(P_LAB(1)**2+P_LAB(2)**2)) IF(P_LAB(2).LT.0.) PHI_P=6.2831853-PHI_P IF ((PHI_P.LT.6.2831853+PHIMIN).AND.(PHI_P.GT.PHIMAX)) GOTO 180 C Neutral pion: I_PARTICLE = 7 CALL UGPKIN(I_PARTICLE,3,BETA,P_LAB,PLB,THETA_PI,GEOMWT) C Find the weight IP = 0 CALL EPC_O (IP,EBEAM,P_LAB(5),THETA_PI,SIGMA_PI) C ns EL_WEIGHT = SIGMA_PI/TOT_MOTT_CROSS/GEOMWT EL_WEIGHT = SIGMA_PI WEIGHT_N=EL_WEIGHT IF (LPTRAK(I_PARTICLE)) THEN CALL GSKINE(PLB,I_PARTICLE,NVERT,0,0,NTT) NUBUF=2 UBUF(1)=0 UBUF(2)=1 IADR=0 CALL GSKINU(NTT,NUBUF,UBUF,IADR) IF (NTT.EQ.0) GO TO 220 ENDIF C C Proton: I_PARTICLE = 14 CALL UGPKIN(I_PARTICLE,2,BETA,P_LAB,PLB,THETA_P,GEOMWT) IF (LPTRAK(I_PARTICLE)) THEN CALL GSKINE(PLB,I_PARTICLE,NVERT,0,0,NTT) NUBUF=2 UBUF(1)=0 UBUF(2)=1 IADR=0 CALL GSKINU(NTT,NUBUF,UBUF,IADR) IF (NTT.EQ.0) GO TO 220 ENDIF C C Electron: I_PARTICLE = 3 CALL UGPKIN(I_PARTICLE,1,BETA,P_LAB,PLB,THETA_E,GEOMWT) IF (LPTRAK(I_PARTICLE)) THEN CALL GSKINE(PLB,I_PARTICLE,NVERT,0,0,NTT) NUBUF=2 UBUF(1)=0 UBUF(2)=1 IADR=0 CALL GSKINU(NTT,NUBUF,UBUF,IADR) IF (NTT.EQ.0) GO TO 220 ENDIF C RETURN C C ****************************************************** C * This part for SHORT RANGE CORRELATION * C * in deuteron * C * e+d ==> e+p+n * C ****************************************************** C 601 continue C Calculate and store vertex CALL GSVERT (VTX,0,0,0,0,NVERT) IF (NVERT.EQ.0) GO TO 230 ** ** Random number precalculation: ** IF(MOD(I_EVENT,10000).EQ.0) PRINT *,I_EVENT IF(I_EVENT.EQ.0) THEN XLOW=300. XHIGH=1000. CALL FUNLXP(FERMI_DIST,FSPACE,XLOW,XHIGH) ENDIF ** ** Generate Fermi momenta: ** I_EVENT = I_EVENT+1 c CALL FUNLUX(FSPACE,PT1,1) c PT1=PT1/1000. c PRINT *, PT1 611 PT1=0.3+0.7*RANDOM() CALL RAN3D(PT1_X,PT1_Y,PT1_Z,PT1) PR1_PX_LAB=PT1_X PR1_PY_LAB=PT1_Y PR1_PZ_LAB=PT1_Z AN_PX_LAB=-PT1_X AN_PY_LAB=-PT1_Y AN_PZ_LAB=-PT1_Z ** ** Fermi motion of the neutron spectator: ** AN_P_LAB=SQRT(AN_PX_LAB**2+AN_PY_LAB**2+AN_PZ_LAB**2) AN_E_LAB=SQRT(MNEUTRON**2+AN_P_LAB**2) T_NEUTRON=AN_E_LAB-MNEUTRON c** c** Fermi motion of the proton specator: c** PR1_P_LAB=SQRT(PR1_PX_LAB**2+PR1_PY_LAB**2+PR1_PZ_LAB**2) PR1_E_LAB=SQRT(MPROTON**2+PR1_P_LAB**2) THETA_PR1=ACOS(PR1_PZ_LAB/PR1_P_LAB) IF(THETA_PR1.LT.0.) THETA_PR1=3.1415927+THETA_PR1 PHI_PR1=ATAN2(PR1_PX_LAB,PR1_PY_LAB) IF(PHI_PR1.LT.0.) PHI_PR1=6.2831853+PHI_PR1 * * Kinetic energy taken by proton: * T_PROTON=PR1_E_LAB-MPROTON ** ** Some of the constans needed: ** Q_REACTION=T_PROTON+T_NEUTRON C_EN=EBEAM+AN_E_LAB-Q_REACTION ** ** The other constans needed: ** V_EN_X=AN_PX_LAB V_EN_Y=AN_PY_LAB V_EN_Z=EBEAM+AN_PZ_LAB V_EN=SQRT(V_EN_X**2+V_EN_Y**2+V_EN_Z**2) IF(V_EN_Z/V_EN.GT.1.) THEN PRINT *, V_EN_Z,V_EN V_EN_Z=V_EN ENDIF THETA_V_EN = ACOS(V_EN_Z/V_EN) PHI_V_EN = ACOS( V_EN_X/SQRT(V_EN_X**2+V_EN_Y**2)) IF(V_EN_Y.LT.0.) PHI_V_EN=6.2831853-PHI_V_EN V_X=SIN(THETA_V_EN)*COS(PHI_V_EN) V_Y=SIN(THETA_V_EN)*SIN(PHI_V_EN) V_Z=COS(THETA_V_EN) ** ** Scattered electron energy and momentum: ** THETA_EL = 0.0002+3.1413927*RANDOM() PHI_EL=6.2831853*RANDOM() E_X_LAB=SIN(THETA_EL)*COS(PHI_EL) E_Y_LAB=SIN(THETA_EL)*SIN(PHI_EL) E_Z_LAB=COS(THETA_EL) c CALL RAN3D(E_X_LAB,E_Y_LAB,E_Z_LAB,1.) c THETA_EL=ACOS(E_Z_LAB) c PHI_EL=ACOS(E_X_LAB/SQRT(E_X_LAB**2+E_Y_LAB**2)) c IF(E_Y_LAB.LT.0.) PHI_EL=6.2831853-PHI_EL CS_EN=E_X_LAB*V_X+E_Y_LAB*V_Y+E_Z_LAB*V_Z OO2=2.*(V_EN*CS_EN-C_EN) IF(OO2.NE.0.) THEN ELECTRON_E_LAB=(V_EN**2-C_EN**2+MNEUTRON**2)/OO2 ENDIF IF(ELECTRON_E_LAB.LT.0..OR.ELECTRON_E_LAB.GT.EBEAM) GOTO 611 ELECTRON_PX_LAB=ELECTRON_E_LAB*SIN(THETA_EL)*COS(PHI_EL) ELECTRON_PY_LAB=ELECTRON_E_LAB*SIN(THETA_EL)*SIN(PHI_EL) ELECTRON_PZ_LAB=ELECTRON_E_LAB*COS(THETA_EL) ELECTRON_P_LAB =ELECTRON_E_LAB RES1=THETA_EL RES2=PHI_EL RES3=ELECTRON_P_LAB ** ** Scattered neutron energy and momentum: ** AN_PX_LAB=0.-ELECTRON_PX_LAB-PR1_PX_LAB AN_PY_LAB=0.-ELECTRON_PY_LAB-PR1_PY_LAB AN_PZ_LAB=EBEAM-ELECTRON_PZ_LAB-PR1_PZ_LAB AN_P_LAB =SQRT(AN_PX_LAB**2+AN_PY_LAB**2+AN_PZ_LAB**2) AN_E_LAB=SQRT(MNEUTRON**2+AN_P_LAB**2)-MNEUTRON THETA_PRN = ACOS(AN_PZ_LAB/AN_P_LAB) IF(THETA_PRN.LT.0.) THETA_PRN=180.+THETA_PRN PHI_PRN = ATAN2(AN_PY_LAB,AN_PX_LAB) IF(PHI_PRN.LT.0.) PHI_PRN=360.+PHI_PRN * * Mott cross section: * CROSS_MOTT | =XXMOTT(EBEAM,ELECTRON_P_LAB,THETA_EL,MNEUTRON,1.) Q__2=4.*EBEAM*ELECTRON_P_LAB*SIN(THETA_EL/2.)**2 T2THE = TAN(THETA_EL/2.)*TAN(THETA_EL/2.) ! ! Galster or dipole parameterzation : ! G_D_v=1./(1.+Q__2/0.71)**2 ! ! Proton and neutron form factors : ! TAU=Q__2/4./MPROTON**2 FAC= 1./(1.+TAU) G_P_E = G_D_v S_n =1./(1.+5.6*Tau) G_P_M = 2.7928456*G_D_v G_N_M = -1.91304184*G_D_v G_N_EL = - 1.91304184 * Tau * G_D_v * S_n * * Cross section: * SIGMA_DIPOLE_N = CROSS_MOTT* | (G_N_EL**2*FAC+TAU*G_N_M**2*(FAC+2.*T2THE)) WEIGHT_N=SIGMA_DIPOLE_N ** Check on PROTON1: THETA_P = ACOS(PT1_Z/PT1) IF(THETA_P.LT.0.) THETA_P=3.1415927+THETA_P IF((THETA_P.LT.THETAMIN).OR.(THETA_P.GT.THETAMAX)) GOTO 611 PHI_P = ACOS(PT1_X/SQRT(PT1_X**2+PT1_Y**2)) IF(PT1_Y.LT.0.) PHI_P=6.2831853-PHI_P C IF((PHI_P.LT.6.2831853+PHIMIN).AND.(PHI_P.GT.PHIMAX)) GOTO 611 ** Follow the Proton: I_PARTICLE = 14 PLB(1)=PT1_X PLB(2)=PT1_Y PLB(3)=PT1_Z P_LAB(1)=PLB(1) P_LAB(2)=PLB(2) P_LAB(3)=PLB(3) P_LAB(5)=SQRT(PLB(1)**2+PLB(2)**2+PLB(3)**2) P_LAB(4)=SQRT(P_LAB(5)**2+MPROTON**2) IF(LPTRAK(I_PARTICLE)) THEN CALL GSKINE(PLB,I_PARTICLE,NVERT,0,0,NTT) NUBUF=2 UBUF(1)=0 UBUF(2)=1 IADR=0 CALL GSKINU(NTT,NUBUF,UBUF,IADR) IF (NTT.EQ.0) GO TO 220 ENDIF ** Follow the Neutron: I_PARTICLE = 13 P_LAB(1)=AN_PX_LAB P_LAB(2)=AN_PY_LAB P_LAB(3)=AN_PZ_LAB PLB(1)=P_LAB(1) PLB(2)=P_LAB(2) PLB(3)=P_LAB(3) P_LAB(5)=SQRT(PLB(1)**2+PLB(2)**2+PLB(3)**2) P_LAB(4)=SQRT(P_LAB(5)**2+MNEUTRON**2) IF(LPTRAK(I_PARTICLE)) THEN CALL GSKINE(PLB,I_PARTICLE,NVERT,0,0,NTT) NUBUF=2 UBUF(1)=0 UBUF(2)=1 IADR=0 CALL GSKINU(NTT,NUBUF,UBUF,IADR) IF (NTT.EQ.0) GO TO 220 ENDIF ** Follow the Electron: I_PARTICLE = 3 CALL UGPKIN(I_PARTICLE,1,BETA,P_LAB,PLB,THETA_E,GEOMWT) IF(LPTRAK(I_PARTICLE)) THEN CALL GSKINE(PLB,I_PARTICLE,NVERT,0,0,NTT) NUBUF=2 UBUF(1)=0 UBUF(2)=1 IADR=0 CALL GSKINU(NTT,NUBUF,UBUF,IADR) IF (NTT.EQ.0) GO TO 220 ENDIF RETURN C C ****************************************************** C * This part for SHORT RANGE CORRELATION * C * in deuteron * C * e+3He ==> e+p+p+n * C ****************************************************** 701 continue C Calculate and store vertex CALL GSVERT (VTX,0,0,0,0,NVERT) IF (NVERT.EQ.0) GO TO 230 ** ** Random number precalculation: ** IF(I_EVENT.EQ.0) THEN XLOW=300. XHIGH=1000. CALL FUNLXP(FERMI_DIST,FSPACE,XLOW,XHIGH) ENDIF ** ** Generate Fermi momenta: ** I_EVENT = I_EVENT+1 IF(MOD(I_EVENT,1000).EQ.0) PRINT *, I_EVENT CCC CALL FUNLUX(FSPACE,PT1,1) !to generate L numbers XRAN CCC PT1=PT1/1000. 711 PT1=0.30+0.6*RANDOM() PT2=0.30+0.6*RANDOM() CALL RAN3D(PT1_X,PT1_Y,PT1_Z,PT1) CALL RAN3D(PT2_X,PT2_Y,PT2_Z,PT2) PR1_PX_LAB=PT1_X PR1_PY_LAB=PT1_Y PR1_PZ_LAB=PT1_Z PR2_PX_LAB=PT2_X PR2_PY_LAB=PT2_Y PR2_PZ_LAB=PT2_Z AN_PX_LAB=-PT1_X-PT2_X AN_PY_LAB=-PT1_Y-PT2_Y AN_PZ_LAB=-PT1_Z-PT2_Z ** ** Fermi motion of the neutron spectator: ** AN_P_LAB =SQRT(AN_PX_LAB**2+AN_PY_LAB**2+AN_PZ_LAB**2) AN_E_LAB=SQRT(MNEUTRON**2+AN_P_LAB**2) T_NEUTRON=AN_E_LAB-MNEUTRON c** c** Fermi motion of the proton 1 specator: c** PR1_P_LAB=SQRT(PR1_PX_LAB**2+PR1_PY_LAB**2+PR1_PZ_LAB**2) PR1_E_LAB=SQRT(MPROTON**2+PR1_P_LAB**2) THETA_PR1 = ACOS(PR1_PZ_LAB/PR1_P_LAB) IF(THETA_PR1.LT.0.) THETA_PR1=3.1415927+THETA_PR1 PHI_PR1 = ATAN2(PR1_PX_LAB,PR1_PY_LAB) IF(PHI_PR1.LT.0.) PHI_PR1=6.2831853+PHI_PR1 c** c** Fermi motion of the proton 2 specator: c** PR2_P_LAB=SQRT(PR2_PX_LAB**2+PR2_PY_LAB**2+PR2_PZ_LAB**2) PR2_E_LAB=SQRT(MPROTON**2+PR2_P_LAB**2) THETA_PR2 = ACOS(PR2_PZ_LAB/PR2_P_LAB) IF(THETA_PR2.LT.0.) THETA_PR2=3.1415927+THETA_PR2 PHI_PR2 = ATAN2(PR2_PX_LAB,PR2_PY_LAB) IF(PHI_PR2.LT.0.) PHI_PR2=6.2831853+PHI_PR2 * * Kinetic energy taken by two protons: * T_PROTON1=PR1_E_LAB-MPROTON T_PROTON2=PR2_E_LAB-MPROTON T_TWO_PROTONS=T_PROTON1+T_PROTON2 ** ** One of the constans needed: ** Q_REACTION=T_TWO_PROTONS+T_NEUTRON C_EN=EBEAM+AN_E_LAB-Q_REACTION ** ** The other constans needed: ** V_EN_X=AN_PX_LAB V_EN_Y=AN_PY_LAB V_EN_Z=EBEAM+AN_PZ_LAB V_EN=SQRT(V_EN_X**2+V_EN_Y**2+V_EN_Z**2) IF(V_EN_Z/V_EN.GT.1.) THEN PRINT *, V_EN_Z,V_EN V_EN_Z=V_EN ENDIF THETA_V_EN = ACOS(V_EN_Z/V_EN) PHI_V_EN = ACOS( V_EN_X/SQRT(V_EN_X**2+V_EN_Y**2)) IF(V_EN_Y.LT.0.) PHI_V_EN=6.2831853-PHI_V_EN V_X=SIN(THETA_V_EN)*COS(PHI_V_EN) V_Y=SIN(THETA_V_EN)*SIN(PHI_V_EN) V_Z=COS(THETA_V_EN) ** ** Scattered electron energy and momentum: ** THETA_EL = 0.0002+3.1413927*RANDOM() PHI_EL=6.2831853*RANDOM() E_X_LAB=SIN(THETA_EL)*COS(PHI_EL) E_Y_LAB=SIN(THETA_EL)*SIN(PHI_EL) E_Z_LAB=COS(THETA_EL) CS_EN=E_X_LAB*V_X+E_Y_LAB*V_Y+E_Z_LAB*V_Z OO2=2.*(V_EN*CS_EN-C_EN) IF(OO2.NE.0.) THEN ELECTRON_E_LAB=(V_EN**2-C_EN**2+MNEUTRON**2)/OO2 ENDIF IF(ELECTRON_E_LAB.LT.0..OR.ELECTRON_E_LAB.GT.EBEAM) GOTO 711 ELECTRON_PX_LAB=ELECTRON_E_LAB*SIN(THETA_EL)*COS(PHI_EL) ELECTRON_PY_LAB=ELECTRON_E_LAB*SIN(THETA_EL)*SIN(PHI_EL) ELECTRON_PZ_LAB=ELECTRON_E_LAB*COS(THETA_EL) ELECTRON_P_LAB =ELECTRON_E_LAB RES1=THETA_EL RES2=PHI_EL RES3=ELECTRON_P_LAB ** ** Scattered neutron energy and momentum: ** AN_PX_LAB=0.-ELECTRON_PX_LAB-PR2_PX_LAB-PR1_PX_LAB AN_PY_LAB=0.-ELECTRON_PY_LAB-PR2_PY_LAB-PR1_PY_LAB AN_PZ_LAB=EBEAM-ELECTRON_PZ_LAB-PR2_PZ_LAB-PR1_PZ_LAB AN_P_LAB =SQRT(AN_PX_LAB**2+AN_PY_LAB**2+AN_PZ_LAB**2) AN_E_LAB=SQRT(MNEUTRON**2+AN_P_LAB**2)-MNEUTRON THETA_PRN = ACOS(AN_PZ_LAB/AN_P_LAB) IF(THETA_PRN.LT.0.) THETA_PRN=180.+THETA_PRN PHI_PRN = ATAN2(AN_PY_LAB,AN_PX_LAB) IF(PHI_PRN.LT.0.) PHI_PRN=360.+PHI_PRN * * Mott cross section: * CROSS_MOTT | =XXMOTT(EBEAM,ELECTRON_P_LAB,THETA_EL,MNEUTRON,1.) Q__2=4.*EBEAM*ELECTRON_P_LAB*SIN(THETA_EL/2.)**2 T2THE = TAN(THETA_EL/2.)*TAN(THETA_EL/2.) ! ! Galster or dipole parameterzation : ! G_D_v=1./(1.+Q__2/0.71)**2 ! ! Proton and neutron form factors : ! TAU=Q__2/4./MPROTON**2 FAC= 1./(1.+TAU) G_P_E = G_D_v S_n =1./(1.+5.6*Tau) G_P_M = 2.7928456*G_D_v G_N_M = -1.91304184*G_D_v G_N_EL = - 1.91304184 * Tau * G_D_v * S_n * * Cross section: * SIGMA_DIPOLE_N = CROSS_MOTT* | (G_N_EL**2*FAC+TAU*G_N_M**2*(FAC+2.*T2THE)) WEIGHT_N=SIGMA_DIPOLE_N ** Check on PROTON1 and PROTON2: THETA_P = ACOS(PT1_Z/PT1) IF(THETA_P.LT.0.) THETA_P=3.1415927+THETA_P IF((THETA_P.LT.THETAMIN).OR.(THETA_P.GT.THETAMAX)) GOTO 711 PHI_P = ACOS(PT1_X/SQRT(PT1_X**2+PT1_Y**2)) IF(PT1_Y.LT.0.) PHI_P=6.2831853-PHI_P C IF((PHI_P.LT.6.2831853+PHIMIN).AND.(PHI_P.GT.PHIMAX)) GOTO 711 THETA_P = ACOS(PT2_Z/PT2) IF(THETA_P.LT.0.) THETA_P=3.1415927+THETA_P IF((THETA_P.LT.THETAMIN).OR.(THETA_P.GT.THETAMAX)) GOTO 711 PHI_P = ACOS(PT2_X/SQRT(PT2_X**2+PT2_Y**2)) IF(PT2_Y.LT.0.) PHI_P=6.2831853-PHI_P C IF((PHI_P.LT.6.2831853+PHIMIN).AND.(PHI_P.GT.PHIMAX)) GOTO 711 ** Follow the Proton: I_PARTICLE = 14 PLB(1)=PT1_X PLB(2)=PT1_Y PLB(3)=PT1_Z P_LAB(1)=PLB(1) P_LAB(2)=PLB(2) P_LAB(3)=PLB(3) P_LAB(5)=SQRT(PLB(1)**2+PLB(2)**2+PLB(3)**2) P_LAB(4)=SQRT(P_LAB(5)**2+MPROTON**2) IF(LPTRAK(I_PARTICLE)) THEN CALL GSKINE(PLB,I_PARTICLE,NVERT,0,0,NTT) NUBUF=2 UBUF(1)=0 UBUF(2)=1 IADR=0 CALL GSKINU(NTT,NUBUF,UBUF,IADR) IF (NTT.EQ.0) GO TO 220 ENDIF ** Follow the Proton: I_PARTICLE = 14 PLB(1)=PT2_X PLB(2)=PT2_Y PLB(3)=PT2_Z P_LAB(1)=PLB(1) P_LAB(2)=PLB(2) P_LAB(3)=PLB(3) P_LAB(5)=SQRT(PLB(1)**2+PLB(2)**2+PLB(3)**2) P_LAB(4)=SQRT(P_LAB(5)**2+MPROTON**2) IF(LPTRAK(I_PARTICLE)) THEN CALL GSKINE(PLB,I_PARTICLE,NVERT,0,0,NTT) NUBUF=2 UBUF(1)=0 UBUF(2)=1 IADR=0 CALL GSKINU(NTT,NUBUF,UBUF,IADR) IF (NTT.EQ.0) GO TO 220 ENDIF ** Follow the Neutron: I_PARTICLE = 13 P_LAB(1)=AN_PX_LAB P_LAB(2)=AN_PY_LAB P_LAB(3)=AN_PZ_LAB PLB(1)=P_LAB(1) PLB(2)=P_LAB(2) PLB(3)=P_LAB(3) P_LAB(5)=SQRT(PLB(1)**2+PLB(2)**2+PLB(3)**2) P_LAB(4)=SQRT(P_LAB(5)**2+MNEUTRON**2) IF(LPTRAK(I_PARTICLE)) THEN CALL GSKINE(PLB,I_PARTICLE,NVERT,0,0,NTT) NUBUF=2 UBUF(1)=0 UBUF(2)=1 IADR=0 CALL GSKINU(NTT,NUBUF,UBUF,IADR) IF (NTT.EQ.0) GO TO 220 ENDIF ** Follow the Electron: I_PARTICLE = 3 CALL UGPKIN(I_PARTICLE,1,BETA,P_LAB,PLB,THETA_E,GEOMWT) IF(LPTRAK(I_PARTICLE)) THEN CALL GSKINE(PLB,I_PARTICLE,NVERT,0,0,NTT) NUBUF=2 UBUF(1)=0 UBUF(2)=1 IADR=0 CALL GSKINU(NTT,NUBUF,UBUF,IADR) IF (NTT.EQ.0) GO TO 220 ENDIF RETURN C C ****************************************************** C * This part is for Moller (electron-electron) * C * elastic scattering (-> REACTION 8) * C * Mark Pitt, Oct. 2002; I started with the existing * C * e-p elastic scattering (REACTION 2) and made * C * appropriate modifications for Moller scattering. * C * Kinematic relations and cross section formulae are * C * from Wagner, et al. NIMA294 (1990) 541-548 and * C * Arrington, et al. NIMA311 (1992) 39-48 * C ****************************************************** C 801 continue C Calculate and store vertex CALL GSVERT (VTX,0,0,0,0,NVERT) IF (NVERT.EQ.0) GO TO 230 C PRINT*,'NVERT NTRACK',NVERT,NTRACK C C Call GENBOD to get one event C Input C Number of outgoing particles NP = 2 C Incident Electron mass AMASS(1) = MELEC C Target Electron mass AMASS(2) = MELEC C Weight KGENEV = 1 C Total cm en. TECM = SQRT(AMASS(1)**2+AMASS(2)**2+2.*EBEAM*AMASS(2)) BETA(1) = 0. BETA(2) = 0. BETA(3) = -SQRT(EBEAM**2-AMASS(1)**2)/(EBEAM+AMASS(2)) C If NTT is set to 0., and the first event does not fit the cut, C GSKINE is not called and a "false error" is detected... C NTT = 0. NTT=1. I_EVENT = I_EVENT+1 IF (ISOTROPE.EQ.1) THEN C This part is to have an isotropic distribution in the Lab. C Which is necessary to compute rates... C Each theta is between 0 and PI ; each phi is between 0 C and 2.PI C C Electron E_THETA=RANDOM()*(E_ANMAX-E_ANMIN)+E_ANMIN E_PHI=RANDOM()*(PHIMAX-PHIMIN)+PHIMIN TH_CM=2.*ATAN(TAN(E_THETA)*SQRT((EBEAM+MELEC)/(2.*MELEC))) E_ENERGY=MELEC+(EBEAM-MELEC)*(COS(TH_CM/2.))**2 E_MOMENTUM=SQRT(E_ENERGY**2-MELEC**2) PLB(1)=E_MOMENTUM*COS(E_PHI)*SIN(E_THETA) PLB(2)=E_MOMENTUM*SIN(E_PHI)*SIN(E_THETA) PLB(3)=E_MOMENTUM*COS(E_THETA) Q__2=4*EBEAM*E_MOMENTUM*SIN(E_THETA/2)**2 C IF (LPTRAK(3)) THEN CALL GSKINE(PLB,3,NVERT,0,0,NTT) NUBUF=2 UBUF(1)=0 UBUF(2)=1 IADR=0 CALL GSKINU(NTT,NUBUF,UBUF,IADR) ENDIF C C C calculate Moller cross-section (units are ub/sr) MOL_CROSS=((1+COS(TH_CM))*(3+(COS(TH_CM))**2))**2 MOL_CROSS=MOL_CROSS/(SIN(TH_CM))**4 MOL_CROSS=(1.99E4)*MOL_CROSS WEIGHT_N=MOL_CROSS*SIN(E_THETA) c WEIGHT_N=MOL_CROSS ELSE 31 CALL GENBOD GOOD_EVENT=.FALSE. C The first of the two Moller electrons I_PARTICLE = 3 CALL UGPKIN(I_PARTICLE,1,BETA,P_LAB,PLB,THETA_M1,GEOMWT) IF (P_LAB(1).NE.0) THEN PHI_M1=ATAN(P_LAB(2)/P_LAB(1)) IF (P_LAB(1).LT.0) PHI_M1=PHI_P+3.14 ELSE IF (P_LAB(2).GE.0) PHI_M1=3.14/2 IF (P_LAB(2).LT.0) PHI_M1=-3.14/2 ENDIF IF (LPTRAK(I_PARTICLE)) THEN C use the limits for electrons defined in qweak.bash IF ((THETA_M1.GT.E_ANMIN).AND.(THETA_M1.LT.E_ANMAX).AND. $ (PHI_M1.GT.PHIMIN).AND.(PHI_M1.LT.PHIMAX)) THEN CALL GSKINE(PLB,I_PARTICLE,NVERT,0,0,NTT) IF (NTT.EQ.0) GO TO 220 GOOD_EVENT=.TRUE. NUBUF=2 UBUF(1)=0 UBUF(2)=1 IADR=0 CALL GSKINU(NTT,NUBUF,UBUF,IADR) ENDIF C The above line was changed to (to quiet the compiler): IF (.NOT. GOOD_EVENT) GOTO 31 ENDIF C C The second of the two Moller electrons: I_PARTICLE = 3 C CALL UGPKIN(I_PARTICLE,2,BETA,P_LAB,PLB,THETA_E,GEOMWT) C Find the weight C calculate Moller cross-section (units are ub/sr) TH_CM=2.*ATAN(TAN(THETA_E)*SQRT((EBEAM+MELEC)/(2.*MELEC))) MOL_CROSS=((1+COS(TH_CM))*(3+(COS(TH_CM))**2))**2 MOL_CROSS=MOL_CROSS/(SIN(TH_CM))**4 MOL_CROSS=(1.99E4)*MOL_CROSS WEIGHT_N = MOL_CROSS/WT IF (LPTRAK(I_PARTICLE)) THEN GOOD_EVENT=.FALSE. IF ((THETA_E.GT.E_ANMIN).AND.(THETA_E.LT.E_ANMAX)) THEN CALL GSKINE(PLB,I_PARTICLE,NVERT,0,0,NTT) IF (NTT.EQ.0) GO TO 220 NUBUF=2 UBUF(1)=0 UBUF(2)=1 IADR=0 CALL GSKINU(NTT,NUBUF,UBUF,IADR) GOOD_EVENT=.TRUE. ENDIF ENDIF C Normally, return after elastic event generation. But C in special plot mode, go on to generate an inelastic event C for the plot. If EITHER electron or proton did not satify C the cut, generate a new one IF(.NOT.GOOD_EVENT)THEN C Remove any other particule whose cut was satisfied NTRACK=0 GO TO 31 ELSE Q__2=QQ__2 ENDIF ENDIF RETURN C ************************************* C * This part for inclusive reaction: * C * e + p --> part + ...... * C ************************************* C This reaction is described using Lightbody O'Connell code EPC C User defines which particle is generated. He also defines the C range in energy. Angles are chosen according to cuts in the C data cards. 320 continue C Compute vertex CALL GSVERT (VTX,0,0,0,0,NVERT) IF (NVERT.EQ.0) GO TO 230 C Particle mass I_PARTICLE=I_REACTION-10 CALL GFPART(I_PARTICLE,NAIPART,ITRTYPE,AMASS(1),CHARGE, $ TLIFE,UB,NWB) I_EVENT = I_EVENT+1 IF(MOD(I_EVENT,1000).EQ.0) PRINT *, I_EVENT THETA=RANDOM()*(THETAMAX-THETAMIN)+THETAMIN PHI=RANDOM()*(PHIMAX-PHIMIN)+PHIMIN c If not using the Brems generator, do the following if(itarb.eq.0) then ! not Brems C The minimum energy is particle at rest : AMASS(1) C The upper limit in energy is EBEAM IF(EBEAM.LT.AMASS(1)+CUTHAD) GOTO 220 c P_LAB(5)=RANDOM()*(EBEAM-AMASS(1)-CUTHAD)+AMASS(1)+CUTHAD c P_LAB(4)=SQRT(P_LAB(5)**2-AMASS(1)**2) c The above two lines are incorrect. The epc code generates c cross sections like c d sigma / (dp_lab * d Omega_lab) c NOT c d sigma / (dE_lab * d Omega_lab). c So, this was changed so that the lab momentum is distributed flat: p_min=sqrt((cuthad+amass(1))**2-amass(1)**2) c approximate expression for maximum momentum IF(ITARGET.EQ.3) THEN p_max=sqrt((ebeam-0.01719)**2-amass(1)**2) ELSEIF(ITARGET.EQ.1.AND.I_PARTICLE.EQ.9) THEN p_max=sqrt((ebeam-amass(1))**2-amass(1)**2) ELSE p_max=sqrt(ebeam**2-amass(1)**2) ENDIF p_lab(4)=random()*(p_max-p_min)+p_min p_lab(5)=sqrt(p_lab(4)**2+amass(1)**2) PLB(1)=P_LAB(4)*SIN(THETA)*COS(PHI) PLB(2)=P_LAB(4)*SIN(THETA)*SIN(PHI) PLB(3)=P_LAB(4)*COS(THETA) IF (LPTRAK(I_PARTICLE)) THEN CALL GSKINE(PLB,I_PARTICLE,NVERT,0,0,NTT) IF (NTT.EQ.0) GO TO 220 NUBUF=2 UBUF(1)=0 UBUF(2)=1 IADR=0 CALL GSKINU(NTT,NUBUF,UBUF,IADR) IF (NTT.EQ.0) GO TO 220 Q__2=AMASS(1)**2+MELEC**2-2*P_LAB(4)*EBEAM Q__2=Q__2+2*P_LAB(5)*EBEAM*COS(THETA) Q__2=-Q__2 C Find the weight C IP= 1 ==> proton C IP=-1 ==> neutron C IP= 2 ==> positive pion C IP=-2 ==> negative pion C IP= 0 ==> neutral pion IF (I_PARTICLE.EQ.14) THEN IP=1 ELSEIF (I_PARTICLE.EQ.13) THEN IP=-1 ELSEIF (I_PARTICLE.EQ.8) THEN IP=2 ELSEIF (I_PARTICLE.EQ.9) THEN IP=-2 ELSEIF (I_PARTICLE.EQ.7) THEN IP=0 ELSE GO TO 220 ENDIF if(itarget.eq.1) then ! hydrogen target CALL EPC_O(IP,EBEAM,P_LAB(4),THETA,SIGMA_PI,1,0,0) elseif(itarget.eq.2) then ! deuterium target CALL EPC_O(IP,EBEAM,P_LAB(4),THETA,SIGMA_PI,1,1,mode) elseif(itarget.eq.3) then ! neutron CALL EPC_O(IP,EBEAM,P_LAB(4),THETA,SIGMA_PI,0,1,mode) else c use Jeff's single nucleon pi- cross section instead c (can easily swap in EPC as desired by removing this branch) call delta_jeff(4,ebeam*1000., | (p_lab(5)-amass(1))*1000.,theta,sigma_pi) c convert to dsigma / (domegapi dppi) by multiplying by ppi/Epi Jacobian. sigma_pi=sigma_pi*p_lab(4)/p_lab(5) endif endif ! IPTY else ! Brems generator branch if(LPTRAK(I_PARTICLE)) THEN ! why bother if not tracking? c In this branch, p_nu is actually the lab PHOTON momentum and theta c is the CM PION momentum. c This is necessary due to the form of the cross section. c PHI is the CM or LAB (equivalent) XY plane angle. c We will let GEANT worry about if the pion hasn't enough momentum c to get tracked ... will generate some error messages but who cares. p_min=0.152 ! pion threshold p_max=ebeam-0.001 ! take an MeV off the beam for NAN safety p_nu=random()*(p_max-p_min)+p_min c calculate weight c see maid_g0.F or maid2000.F for more info c note! this is untested for all except pi- branch!!! IF (I_PARTICLE.EQ.14) THEN ! proton ISO=4 ! via delta0 ELSEIF (I_PARTICLE.EQ.13) THEN ! neutron ISO=3 ! via delta+ ELSEIF (I_PARTICLE.EQ.8) THEN ! pi+ ISO=3 ! via delta+ ELSEIF (I_PARTICLE.EQ.9) THEN ! pi- ISO=4 ! via delta0 ELSEIF (I_PARTICLE.EQ.7) THEN ! pi0 ISO=1 ! via delta+ ELSE GO TO 220 ENDIF if(abs(ipty-2).ne.1)then ! vertex not at end of target if(itarget.eq.1.or.itarget.eq.2)then call pi_brems(iso,1,itarget,bxtarg(3)*2., c ebeam*1000.,p_nu*1000., c theta,sigma_pi,ppi,thetapi) else stop "UNKNOWN TARGET TYPE" endif elseif(ipty.eq.1)then ! upstream Al window c For G0, there are the front and back He windows upstream of the target. c They are each 5 mil thick. call pi_brems2(iso,1,2,0.0,0.01905, c ebeam*1000.,p_nu*1000., c theta,sigma_pi,ppi,thetapi) if(iso.eq.3)then sigma_pi=sigma_pi*14. ! normalize per nucleus else sigma_pi=sigma_pi*13. ! normalize per nucleus endif else ! downstream Al window gets slammed by brems call pi_brems2(iso,1,2,bxtarg(3)*2.,0.0254, c ebeam*1000.,p_nu*1000., c theta,sigma_pi,ppi,thetapi) if(iso.eq.3)then sigma_pi=sigma_pi*14. ! normalize per nucleus else sigma_pi=sigma_pi*13. ! normalize per nucleus endif endif c convert back to GeV and make pion four-vector p_lab(4)=ppi/1000. PLB(1)=P_LAB(4)*SIN(THETAPI)*COS(PHI) PLB(2)=P_LAB(4)*SIN(THETAPI)*SIN(PHI) PLB(3)=P_LAB(4)*COS(THETAPI) p_lab(5)=sqrt(p_lab(4)**2+amass(1)**2) c insert pion momentum into event (copied from other branch) CALL GSKINE(PLB,I_PARTICLE,NVERT,0,0,NTT) IF (NTT.EQ.0) GO TO 220 NUBUF=2 UBUF(1)=0 UBUF(2)=1 IADR=0 CALL GSKINU(NTT,NUBUF,UBUF,IADR) IF (NTT.EQ.0) GO TO 220 Q__2=AMASS(1)**2+MELEC**2-2*P_LAB(4)*EBEAM Q__2=Q__2+2*P_LAB(5)*EBEAM*COS(THETAPI) Q__2=-Q__2 q__2=p_nu ! (make q__2=p_nu for tests) endif ! lptrak if endif ! brems if C The weight is proportional to SIGMA_PI WEIGHT_N=SIGMA_PI ! gives weight_n in ub return C Error returned by GSKINE. 220 WRITE (CHMAIL,210) I_REACTION,I_PARTICLE 210 FORMAT (' GSKINE returned an error for reaction',I2, $ ' and particle type',I3,'.') CALL GMAIL(0,0) CALL ERSTOP('Error setting up kinematics in GUKINE',0) C Error returned by GSVERT 230 WRITE (CHMAIL,240) I_REACTION 240 FORMAT (' GSVERT returned an error for reaction',I2,'.') CALL GMAIL(0,0) CALL ERSTOP('Error setting up vertex in GSVERT',0) END c------------------------------------------------------------------------- c** c** Mott scatering : c** FUNCTION XXMOTT(EI,EF,THE,XMT,Z) C Function to calculate Mott cross section in mikrobarns/sr C including recoil. REAL EI,EF,XMT,Z,XSEC,EIM,EFM,XMTM EIM=1000.*EI EFM=1000.*EF XMTM=1000.*XMT CTH=COS(THE/2.) STH=SIN(THE/2.) c ETA=1.+2.*EI/XMTM*STH*STH c XSEC=(Z*1.44/2./EI*CTH/STH/STH)**2/ETA XSEC=(Z*1.44/2./EIM*CTH/STH/STH)**2*(EFM/EIM) XXMOTT=XSEC*1.E4 RETURN END c** c** Fermi distribution for Deuteron: c** FUNCTION FERMI_DIST(X) C Function to calculate Fermi distribution FERMI_DIST=2.*10.**10*(1./(2088.+x**2)/(56460.+x**2))**2 | +10.*10.**10*(x**2/(2088.+x**2)/(96549.+x**2)**2)**2 FERMI_DIST=x**2*FERMI_DIST RETURN END