C PROGRAM EPCV(INPUT,TAPE5=INPUT) C ELECTROPRODUCTION YIELDS OF NUCLEONS AND PIONS C WRITTEN BY J.S. O'CONNELL AND J.W. LIGHTBODY, JR. C NATIONAL BUREAU OF STANDARDS C SEPTEMBER 1987 ! all the variables finishing by al were added to the original code ! and all the comment with ! ! and all the non-capital text. ! CEBAF A. Deur 07/30/98 ! this code compute the counting rates at CEBAF for e94-010 SUBROUTINE WISER(Z,N,PART,E_IN,P_IN,TH_IN,radlen_IN,XSEC) ! Z, N are target info ! PART is the particle whose cross section is to be calcualted ! E_IN is beam energy in GeV ! P_IN is scattered particle momentum in GeV ! TH_IN is scattering angle in rad ! radlen_IN is total radiation length in % ! XSEC is the output cross section in nbarn/GeV/sr, per nucleon ! IMPORTANT: to get the cross section per nuclei, times this by (A)**0.8 ! ! X. Zheng June 2006 IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER Z,N CHARACTER*3 PART COMMON/QD/QDF COMMON/DEL/IP ! Note: mass are used in MeV only in wiser_sub.f and pos_wiser_sub.f DATA PI/3.1416/,AM/938.28/,AMD/1875.63/,AMPI/139.6/ ! c parameter (q3al=1.6,solang3al=0.0039,convbarn=1.e4) ! solang HRS acceptance c parameter (acc3al=8./100.) !HRS momentum acceptance logical WDEBUG c WDEBUG=.true. WDEBUG=.false. E1=E_IN*1000. P=P_IN*1000. THP=TH_IN*180./PI rad_len=radlen_IN IA=Z+N if (WDEBUG) then write(*,'("in WISER:",3I4,4F8.4)') Z,N,IA,E1, P, THP, targL endif IF(IA.EQ.2.OR.IA.EQ.3.OR.IA.EQ.4.OR.IA.EQ.12.OR.IA.EQ.16)THEN ELSEIF(IA.EQ.1)THEN PRINT *, ' NUCLEON CASE' ELSE if (WDEBUG) PRINT 107 107 FORMAT(' NO SPECIFIC SPECTRAL FUNCTION FOR THIS NUCLEUS') if (WDEBUG) PRINT 108 108 FORMAT(' A = 16 SPECTRAL FUNCTION WILL BE USED') ENDIF C 'AN' IS EFFECTIVE NUMBER OF NUCLEONS FOR PION PRODUCTION IF(PART.EQ.'P')THEN AN=N/3.+2.*Z/3. IP=1 ITYPE=5 ELSEIF(PART.EQ.'N')THEN AN=Z/3.+2.*N/3. IP=-1 ELSEIF(PART.EQ.'PI+')THEN AN=Z/3. IP=2 ITYPE=1 ELSEIF(PART.EQ.'PI-')THEN AN=N/3. IP=2 ITYPE=2 ELSEIF(PART.EQ.'PI0')THEN AN=2.*(N+Z)/3. IP=2 ELSE STOP ENDIF if (WDEBUG) > write(*,'("AN=",F4.1," IP=",I2," ITYPE=",I2)') AN,IP,ITYPE IF(IA.EQ.1)THEN if (WDEBUG) PRINT 114 114 FORMAT(' ONLY DELTA MECHANISM INCLUDED') ELSEIF(ABS(IP).EQ.1)THEN IF(IA.GT.1.AND.IA.LT.5)THEN DLF=IA ELSE DLF=7. ENDIF if (WDEBUG) PRINT 102 102 FORMAT(' ENTER LEVINGER FACTOR OR 0 FOR DEFAULT VALUE') READ *,AL IF(AL.EQ.0.)AL=DLF QDF=AL*N*Z/FLOAT(IA) ELSE ENDIF PTP=P IF(ABS(IP).EQ.1)THEN E=SQRT(PTP**2+AMP**2) TP=PTP**2/(E+AMP) AJ=PTP/E D2QD=0. D2QF=0. ELSE E=SQRT(PTP**2+AMPI**2) TP=PTP**2/(E+AMPI) AJ=PTP/E D2QD=0. D2QF=0. ENDIF if (WDEBUG) print *,'E=',E,TP,AJ,D2QD,D2QF D2DEL=0. IF(ABS(IP).EQ.2.OR.IP.EQ.0)THEN if (WDEBUG) print *,'rad_len=',rad_len if (WDEBUG) print *,'calling WISER_ALL_SIG', > E1,PTP,THP,RAD_LEN,ITYPE,TOTAL Call WISER_ALL_SIG(E1,PTP,THP,RAD_LEN,ITYPE,TOTAL) ! TOTAL in nanobarn/GeV*str, for proton target XSEC=TOTAL ! should be timed by "effective number of nucleons per nuclei"??? ! or (N+Z)??? c XSEC=XSEC*((N+Z)**0.8) ! now is for per nuclei if (WDEBUG) print *,'WISER end, XSEC=',XSEC ENDIF 999 continue RETURN END C----------------------------------------------------------------------- C----------------------------------------------------------------------- C Xiaochao Zheng, this subroutine added July 07, 2000 DOUBLE PRECISION FUNCTION QUADMO(FUNCT,PLOWER,PUPPER,EPSLON,NLVL) REAL*8 FUNCT,PLOWER,PUPPER,EPSLON INTEGER NLVL INTEGER LEVEL,MINLVL/3/,MAXLVL/24/,IRETRN(50),I REAL*8 VALINT(50,2), VMX(50), RX(50), FMX(50), FRX(50), 1 FMRX(50), ESTRX(50), EPSX(50) REAL*8 R, FL, FML, FM, FMR, FR, EST, ESTL, ESTR, ESTINT,VL, 1 AREA, ABAREA, VM, COEF, ROMBRG, EPS ! IMPLICIT REAL*8 (A-H,O-Z) LEVEL = 0 NLVL = 0 ABAREA = 0.0 VL = PLOWER R = PUPPER FL = FUNCT(VL) FM = FUNCT(0.5*(VL+R)) FR = FUNCT(R) EST = 0.0 EPS = EPSLON 100 LEVEL = LEVEL+1 VM = 0.5*(VL+R) COEF = R-VL IF(COEF.NE.0) GO TO 150 ROMBRG = EST GO TO 300 150 FML = FUNCT(0.5*(VL+VM)) FMR = FUNCT(0.5*(VM+R)) ESTL = (FL+4.0*FML+FM)*COEF ESTR = (FM+4.0*FMR+FR)*COEF ESTINT = ESTL+ESTR AREA=DABS(ESTL)+DABS(ESTR) ABAREA=AREA+ABAREA-DABS(EST) IF(LEVEL.NE.MAXLVL) GO TO 200 NLVL = NLVL+1 ROMBRG = ESTINT GO TO 300 200 IF((DABS(EST-ESTINT).GT.(EPS*ABAREA)).OR. 1 (LEVEL.LT.MINLVL)) GO TO 400 ROMBRG = (1.6D1*ESTINT-EST)/15.0D0 300 LEVEL = LEVEL-1 I = IRETRN(LEVEL) VALINT(LEVEL, I) = ROMBRG GO TO (500, 600), I 400 IRETRN(LEVEL) = 1 VMX(LEVEL) = VM RX(LEVEL) = R FMX(LEVEL) = FM FMRX(LEVEL) = FMR FRX(LEVEL) = FR ESTRX(LEVEL) = ESTR EPSX(LEVEL) = EPS EPS = EPS/1.4 R = VM FR = FM FM = FML EST = ESTL GO TO 100 500 IRETRN(LEVEL) = 2 VL = VMX(LEVEL) R = RX(LEVEL) FL = FMX(LEVEL) FM = FMRX(LEVEL) FR = FRX(LEVEL) EST = ESTRX(LEVEL) EPS = EPSX(LEVEL) GO TO 100 600 ROMBRG = VALINT(LEVEL,1)+VALINT(LEVEL,2) IF(LEVEL.GT.1) GO TO 300 QUADMO = ROMBRG /12.0D0 RETURN END Subroutine WISER_ALL_SIG(E0MM,PMM,THETA_DEG,RAD_LEN,TYPE,SIGMA) !------------------------------------------------------------------------------ ! Calculate pi,K,p cross section for electron beam on a proton target ! IntegrateQs over function WISER_FIT using integration routine QUADMO ! E0 is electron beam energy, OR max of Brem spectra ! P,E is scattered particle momentum,energy ! THETA_DEG is kaon angle in degrees ! RAD_LEN (%)is the radiation length of target, including internal ! (typically 5%) ! = .5 *(target radiation length in %) +5. ! *** =100. IF BREMSTRULUNG PHOTON BEAM OF 1 EQUIVIVENT QUANTA *** ! TYPE: 1 for pi+; 2 for pi-, 3=k+, 4=k-, 5=p, 6=p-bar ! SIGMA is output cross section in nanobars/GeV-str !------------------------------------------------------------------------------ IMPLICIT NONE DOUBLE PRECISION E0MM,PMM,E0,P,THETA_DEG,RAD_LEN,SIGMA INTEGER TYPE COMMON/WISER_ALL/ E,P_COM,COST,P_T,TYPE_COM,PARTICLE,M_X,U_MAN REAL*8 E,P_COM,COST,P_T,M_X,U_MAN C DOUBLE PRECISION E1,P_COM1,COST1,P_T1,M_X1,U_MAN1 INTEGER TYPE_COM,PARTICLE ! Wiser's fit pi+ pi- k+ k- p+ p- REAL*8 A5(6)/-5.49, -5.23, -5.91, -4.45, -6.77, -6.53/ REAL*8 A6(6)/-1.73, -1.82, -1.74, -3.23, 1.90, -2.45/ REAL*8 MASS2(3)/.019488, .2437, .8804/ REAL*8 MASS(3)/.1396, .4973, .9383/ REAL*8 MP/.9383/, MP2/.8804/, RADDEG/.0174533/ REAL*8 M_L,SIG_E REAL*8 E_GAMMA_MIN,WISER_ALL_FIT,QUADMO,E08,EPSILON/.003/ EXTERNAL WISER_ALL_FIT INTEGER N,CHARGE SIG_E = QUADMO(WISER_ALL_FIT,E_GAMMA_MIN,E08,EPSILON,N) * > EXP(A5(TYPE) *M_L) *EXP(A6(TYPE) *P_T**2/E) SIG_E = QUADMO(WISER_ALL_FIT,E_GAMMA_MIN,E08,EPSILON,N) c print *,'in WISER_ALL_SIG',E0MM,PMM,E0,P,THETA_DEG,RAD_LEN,SIGMA P=PMM/1000. E0=E0MM/1000. P_COM = P TYPE_COM = TYPE PARTICLE = (TYPE+1)/2 ! 1= pi, 2= K, 3 =P CHARGE = TYPE -2*PARTICLE +2 ! 1 for + charge, 2 for - charge E08 =E0 E =SQRT(MASS2(PARTICLE) + P**2) COST = COS(RADDEG * THETA_DEG) P_T = P * SIN(RADDEG * THETA_DEG) IF(TYPE.LE.4) THEN !mesons IF(CHARGE.EQ.1) THEN ! K+ n final state M_X = MP ELSE ! K- K+ P final state M_X = MP+ MASS(PARTICLE) ENDIF ELSE ! baryons IF(CHARGE.EQ.1) THEN ! pi p final state M_X = MASS(1) ! pion mass ELSE ! P P-bar P final state M_X = 2.*MP ENDIF ENDIF E_GAMMA_MIN = (M_X**2 -MASS2(PARTICLE ) -MP2+2.*MP*E)/ > (2.*(MP -E +P*COST)) ! WRITE(10,'(''E_GAMMA_MIN='',F10.2,'' p_t='',F8.2)') ! > E_GAMMA_MIN,P_T ! E_GAMMA_MIN = MP *(E + MASS(PARTILCE))/(MP -P*(1.-COST)) * print *,E_GAMMA_MIN IF(E_GAMMA_MIN.GT..1) THEN !Kinematically allowed? M_L = SQRT(P_T**2 + MASS2(PARTICLE)) IF(TYPE.NE.5) THEN ! everything but proton c print *,'calling quadmo:',E_GAMMA_MIN,E08,EPSILON,N c print *,'quadmo=', c > QUADMO(WISER_ALL_FIT,E_GAMMA_MIN,E08,EPSILON,N), c > EXP(A5(TYPE) *M_L), EXP(A6(TYPE) *P_T**2/E) SIG_E = QUADMO(WISER_ALL_FIT,E_GAMMA_MIN,E08,EPSILON,N) * > EXP(A5(TYPE) *M_L) *EXP(A6(TYPE) *P_T**2/E) ELSE ! proton U_MAN = ABS(MP2 + MASS2(PARTICLE) -2.*MP*E) c print *,'quadmo=', c > QUADMO(WISER_ALL_FIT,E_GAMMA_MIN,E08,EPSILON,N), c > EXP(A5(TYPE) *M_L) SIG_E = QUADMO(WISER_ALL_FIT,E_GAMMA_MIN,E08,EPSILON,N) * > EXP(A5(TYPE) *M_L) ENDIF SIGMA = P**2/E * 1000. * RAD_LEN/100. *SIG_E ELSE ! Kinematically forbidden SIGMA = 0. ENDIF c print *,'SIGMA=',SIGMA,SIG_E,RAD_LEN,P,E RETURN END REAL*8 FUNCTION WISER_ALL_FIT(E_GAMMA) !--------------------------------------------------------- ! Calculates pi, k, p cross section for gamma + p -> k ! It is already divided by E_GAMMA, the bremstrulung spectra ! David Wiser's fit from Thesis, eq. IV-A-2 and Table III. ! Can be called from WISER_SIG using integration routine QUADMO ! E,P are KAON energy and momentum ! P_t is KAON transverse momentum ! P_CM is KAON center of mass momentum ! P_CM_L is KAON center of mass longitudinal momentum ! TYPE: 1 for pi+; 2 for pi-, 3=k+, 4=k-, 5=p, 6=p-bar ! E_GAMMA is photon energy. ! Steve Rock 2/21/96 !--------------------------------------------------------- IMPLICIT NONE COMMON/WISER_ALL/ E,P,COST,P_T,TYPE,PARTICLE,M_X,U_MAN REAL*8 E,P,COST,P_T,M_X,U_MAN INTEGER TYPE ! 1 for pi+; 2 for pi-, 3=k+, 4=k-, 5=p, 6=p-bar INTEGER PARTICLE ! 1= pi, 2= K, 3 =P ! Wiser's fit pi+ pi- k+ k- p+ p- REAL*8 A1(6)/566., 486., 368., 18.2, 1.33E5, 1.63E3 / REAL*8 A2(6)/829., 115., 1.91, 307., 5.69E4, -4.30E3 / REAL*8 A3(6)/1.79, 1.77, 1.91, 0.98, 1.41, 1.79 / REAL*8 A4(6)/2.10, 2.18, 1.15, 1.83, .72, 2.24 / REAL*8 A6/1.90/,A7/-.0117/ !proton only REAL*8 MASS2(3)/.019488, .2437, .8804/ REAL*8 MASS(3)/.1396, .4973, .9383/ REAL*8 MP2/.8804/,MP/.9383/, RADDEG/.0174533/ REAL*8 X_R,S,B_CM, GAM_CM, P_CM REAL*8 P_CM_MAX, P_CM_L REAL*8 E_GAMMA !Mandlestam variables S = MP2 + 2.* E_GAMMA * MP !Go to Center of Mass to get X_R B_CM = E_GAMMA/(E_GAMMA+MP) GAM_CM = 1./SQRT(1.-B_CM**2) P_CM_L = -GAM_CM *B_CM *E + > GAM_CM * P * COST P_CM = SQRT(P_CM_L**2 + P_T**2) P_CM_MAX =SQRT (S +(M_X**2-MASS2(PARTICLE))**2/S > -2.*(M_X**2 +MASS2(PARTICLE)) )/2. X_R = P_CM/P_CM_MAX IF(X_R.GT.1.) THEN ! Out of kinematic range WISER_ALL_FIT = 0. ELSEIF(TYPE.NE.5) THEN ! not the proton WISER_ALL_FIT = (A1(TYPE) + A2(TYPE)/SQRT(S)) * > (1. -X_R + A3(TYPE)**2/S)**A4(TYPE)/E_GAMMA ELSE ! special formula for proton WISER_ALL_FIT = ( (A1(TYPE) + A2(TYPE)/SQRT(S)) * > (1. -X_R + A3(TYPE)**2/S)**A4(TYPE) / > (1.+U_MAN)**(A6+A7*S) )/E_GAMMA ENDIF RETURN END