-
- SUBROUTINE fsiini
+
+ SUBROUTINE FSIINI
C---Note:
C-- ICH= 0 (1) if the Coulomb interaction is absent (present);
C-- ISPIN= JJ= 1,2,..,MSPIN denote increasing values of the pair
C-- FD= FDH(LL,JJ), RD= RDH(LL,JJ)= scattering length and effective
C-- radius for each value of the total pair spin S, JJ= 1,..,MSPIN; ;
C-- the corresponding square well parameters EB= EBH(LL,JJ), RB=
-C-- RBH(LL,JJ) (required if NS=1) may be calculated by sear.f;
+C-- RBH(LL,JJ) (required if NS=1) may be calculated by SEAR;
C-- if the effective range approximation is not valid (as is the case,
C-- e.g., for two-pion system) a code for calculation of the scattering
C-- amplitude should be supplemented;
COMMON/FSI_FD/FD(10),RD(10)
COMMON/FSI_C/C(10),AM,AMS,DM
COMMON/FSI_CONS/PI,PI2,SPI,DR,W
- COMPLEX*16 C
+ COMPLEX*16 C
COMMON/FSI_AA/AA
COMMON/FSI_AAPI/AAPI(20,2)/FSI_AAND/AAND(20,4)
- COMMON/FSI_AAPIN/AAPIN(20,2)
COMMON/FSI_SW/RB(10),EB(10),BK(10),CDK(10),SDK(10),
1 SBKRB(10),SDKK(10)
- COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS
-
-c Include 'common_fsi_poc.inc'
-c Include 'common_fsi_prf.inc'
-c Include 'common_fsi_spin.inc'
-c Include 'common_fsi_ach.inc'
-c Include 'common_fsi_ns.inc'
-c Include 'common_fsi_fd.inc'
-c-mlv Include 'common_fsi_c.inc'
-c Include 'common_fsi_cons.inc'
-c Include 'common_fsi_aa.inc'
-c Include 'common_fsi_aapi.inc'
-c Include 'common_fsi_aapin.inc'
-c Include 'common_fsi_sw.inc'
-c Include 'common_fsi_aand.inc'
+ COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS
DIMENSION FDH(30,10),RDH(30,10),EBH(30,10),RBH(30,10)
DIMENSION RHOH(30,10)
INTEGER*4 IERR
C
C--- mass of the first and second particle
- DATA AM1H/.93956563D0,.93827231D0,.93956563D0,3.72737978D0,
- C .13957D0,.13498D0,.13957D0, .93956563D0, .93827231D0,
+ DATA AM1H/.9395656D0,.9382723D0,.9395656D0,3.7294D0,.13957D0,
+ C .13498D0,.13957D0, .9395656D0, .9382723D0,
C 4*.13957D0,4*.493677D0,
- C 2*1.87561339D0,2*2.80892165D0,2*.497672D0,
- C 1.87561339D0,3*.93827231D0,.93956563D0, 2*0.D0/
- DATA AM2H/.93956563D0,.93827231D0,.93827231D0,3.72737978D0,
- C .13957D0,.13498D0,.13957D0, 2*1.87561339D0,
- C 2*.493677D0,2*.93827231D0,
- C 2*.493677D0,2*.93827231D0,
- C 1.87561339D0,3.72737978D0,2.80892165D0,3.72737978D0,
- C 2*.497672D0,2*2.80892165D0,3.72737978D0,
- C 2*1.115684D0,2*0.D0/
-c--------|---------|---------|---------|---------|---------|---------|----------
+ C 2*1.875613D0,2*2.808D0,2*.497672D0,
+ C 1.875613D0,2*.9382723D0, 4*0.D0/
+ DATA AM2H/.9395656D0,.9382723D0,.9382723D0,3.7294D0,.13957D0,
+ C .13498D0,.13957D0, 2*1.875613D0,
+ C 2*.493677D0,2*.9382723D0,
+ C 2*.493677D0,2*.9382723D0,
+ C 1.875613D0,3.7294D0,2.808D0,3.7294D0,
+ C 2*.497672D0,2*2.808D0,3.7294D0, 4*0.D0/
C--- charge of the first and second particle
DATA C1H/0.D0,1.D0,0.D0,2.D0, 1.D0,0.D0,1.D0,0.D0,1.D0,
C 3*1.D0,-1.D0,3*1.D0,-1.D0,
- C 4*1.D0,2*0.D0,4*1.D0,0.D0, 2*0.D0/
+ C 4*1.D0,2*0.D0,3*1.D0, 4*0.D0/
DATA C2H/0.D0,1.D0,1.D0,2.D0,-1.D0,0.D0,3*1.D0,
C -1.D0,3*1.D0,-1.D0,3*1.D0,
- C 1.D0,2.D0,1.D0,2.D0,2*0.D0,2*1.D0,2.D0,2*0.D0,2*0.D0/
+ C 1.D0,2.D0,1.D0,2.D0,2*0.D0,2*1.D0,2.D0, 4*0.D0/
C---MSPIN vs (LL)
- DATA MSPINH/3*2,4*1,2*2,8*1,3,1,2,1,2*1,2*2,1,2*2, 2*0/
+ DATA MSPINH/3*2,4*1,2*2,8*1,3,1,2,1,2*1,2*2,1, 4*0/
C---Spin factors RHO vs (LL,ISPIN)
DATA RHOH/3*.25D0, 4*1.D0, 2*.3333D0, 8*1.D0,
- 1 .1111D0,1.D0,.25D0,1.D0,2*1.D0,
- 1 .3333D0,.25D0,1.D0,2*.25D0, 2*0.D0,
- 2 3*.75D0, 4*0.D0, 2*.6667D0, 8*0.D0,
- 2 .3333D0,.0D0,.75D0,.0D0,2*0.D0,
- 2 .6667D0,.75D0,0.D0,2*.75D0, 2*0.D0,
- 3 17*.0D0,.5556D0,3*0.D0, 7*0.D0,2*0.D0,210*0.D0/
+ C .1111D0,1.D0,.25D0,1.D0,2*1.D0,
+ 1 .3333D0,.25D0,1.D0, 4*0.D0,
+ C 3*.75D0, 4*0.D0, 2*.6667D0, 8*0.D0,
+ C .3333D0,.0D0,.75D0,.0D0,2*0.D0,
+ 2 .6667D0,.75D0,0.D0, 4*0.D0,
+ C 17*.0D0,.5556D0,3*0.D0, 5*0.D0,4*0.D0,210*0.D0/
C---Scattering length FD and effective radius RD in fm vs (LL,ISPIN)
DATA FDH/17.0D0,7.8D0,23.7D0,2230.1218D0,.225D0,.081D0,-.063D0,
- 1 -.65D0,-2.73D0,
- 1 .137D0,-.071D0,-.148D0,.112D0,2*1.D-6,-.360D0,
- 1 2*1.D-6,1.344D0,6*1.D-6,-5.628D0,2.18D0,2.40D0, 2*0.D0,
-cc 2 -10.8D0,2*-5.4D0,4*0.D0,-6.35D0,-11.88D0,8*0.D0,9*0.D0,
- 2 3*-5.4D0,4*0.D0,-6.35D0,-11.88D0,8*0.D0,9*0.D0,
- 2 1.93D0,1.84D0,2*0.D0,
- 3 240*0.D0/
-c--------|---------|---------|---------|---------|---------|---------|----------
+ C -.65D0,-2.73D0,
+ C .137D0,-.071D0,-.148D0,.112D0,2*1.D-6,-.360D0,
+ 1 2*1.D-6,1.344D0,6*1.D-6,-5.628D0, 4*0.D0,
+ C -10.8D0,2*-5.4D0,4*0.D0,-6.35D0,-11.88D0,8*0.D0,9*0.D0,4*0.D0,
+ C 240*0.D0/
DATA RDH/2.7D0,2.8D0,2.7D0,1.12139906D0,-44.36D0,64.0D0,784.9D0,
- 1 477.9D0, 2.27D0, 9*0.D0,-69.973D0, 6*0.D0,3.529D0,
- 1 3.19D0,3.15D0, 2*0.D0,
- 2 3*1.7D0,4*0.D0,2.0D0,2.63D0, 17*0.D0,3.35D0,3.37D0, 2*0.D0,
- 3 240*0.D0/
+c--------|---------|---------|---------|---------|---------|---------|----------
+ C 477.9D0, 2.27D0, 9*0.D0,-69.973D0, 6*0.D0,3.529D0, 4*0.D0,
+ C 3*1.7D0,4*0.D0,2.0D0,2.63D0, 17*0.D0, 4*0.D0, 240*0.D0/
C---Corresponding square well parameters RB (width in fm) and
C-- EB =SQRT(-AM*U) (in GeV/c); U is the well height
DATA RBH/2.545739D0, 2.779789D0, 2.585795D0, 5.023544D0,
- 1 .124673D0, .3925180D0,.09D0, 2.D0, 4.058058D0, 17*0.D0,
- 1 2.252623D0, 2.278575D0, 2*0.D0,
- 2 3*2.003144D0,
- 2 4*0.D0, 2.D0, 4.132163D0, 17*0.D0,
- 2 2.272703D0, 2.256355D0, 2*0.D0,
- 3 240*0.D0/
+ C .124673D0, .3925180D0,.09D0, 2.D0, 4.058058D0, 17*0.D0, 4*0.D0,
+ C 3*2.003144D0,
+ C 4*0.D0, 2.D0, 4.132163D0, 17*0.D0, 4*0.D0, 240*0.D0/
DATA EBH/.1149517D0, .1046257D0, .1148757D0, .1186010D0,
- 1 .7947389D0,2.281208D0,8.7D0,.4D0,.1561219D0,17*0.D0,
- 1 .1013293D0, .1020966D0, 2*0.D0,
- 2 3*.1847221D0,
- 2 4*0.D0, .4D0, .1150687D0, 17*0.D0,
- 2 .09736083D0, .09708310D0, 2*0.D0,
- 3 240*0.D0/
+ C .7947389D0,2.281208D0,8.7D0,.4D0,.1561219D0,17*0.D0,4*0.D0,
+ C 3*.1847221D0,
+ C 4*0.D0, .4D0, .1150687D0, 17*0.D0, 4*0.D0, 240*0.D0/
C=======< constants >========================
W=1/.1973D0 ! from fm to 1/GeV
PI=4*DATAN(1.D0)
PI2=2*PI
SPI=DSQRT(PI)
DR=180.D0/PI ! from radian to degree
-
-c WRITE(*,*)'from C++ to fortran W PI PI2 SPI DR',W,PI,PI2,SPI,DR
-
AC1=1.D10
AC2=1.D10
-C=======< condition de calculs >=============
- NUNIT=11 ! for IBM or HP
-C NUNIT=4 ! for SUN in Prague
-c-mlv CALL readint4(NUNIT,'ITEST ',ITEST)
-c-mlv CALL readint4(NUNIT,'LL ',LL) ! Two-particle system
-c-mlv CALL readint4(NUNIT,'NS ',NS)
-c CALL READ_FILE(NUNIT,'ITEST ',CHAR,ITEST,REAL8,IERR)
-c CALL READ_FILE(NUNIT,'LL ',CHAR,LL,REAL8,IERR)
-c CALL READ_FILE(NUNIT,'NS ',CHAR,NS,REAL8,IERR)
-
-
C---setting particle masses and charges
AM1=AM1H(LL)
C put ICH=0 and substitute the strong amplitude parameters by
C the ones not affected by Coulomb interaction
- IF(ITEST.EQ.0)THEN
+ IF(ITEST.EQ.0)THEN
+
ICH=0
IF(C1*C2.NE.0.D0) ICH=1
IQS=0
IF(C1+AM1.EQ.C2+AM2) IQS=1
I3S=0 ! only this option is available
ISI=1
- I3C=0
- IF(CN*ICH.NE.0.D0) I3C=1
+ I3C=1
+
ENDIF
-
-c IF(ITEST.EQ.1)THEN
-c NS=4 ! SPHER. WAVE
-c ICH=0
-c IQS=1
-c ISI=0
-c I3C=0
-
-
-c-mlv CALL readint4(NUNIT,'ICH ',ICH)
-c-mlv CALL readint4(NUNIT,'IQS ',IQS)
-c-mlv CALL readint4(NUNIT,'ISI ',ISI)
-c-mlv CALL readint4(NUNIT,'I3C ',I3C)
-c CALL READ_FILE(NUNIT,'ICH ',CHAR,ICH,REAL8,IERR)
-c CALL READ_FILE(NUNIT,'IQS ',CHAR,IQS,REAL8,IERR)
-c CALL READ_FILE(NUNIT,'ISI ',CHAR,ISI,REAL8,IERR)
-c CALL READ_FILE(NUNIT,'I3C ',CHAR,I3C,REAL8,IERR)
-c ENDIF
-
- write(*,*)'====itest ll ich iqs isi i3c===',itest,ll,ich, iqs, isi, i3c
-
+c23456
+ write(*,*)'FSIINI ITEST ich iqs i3s isi i3c',ITEST,
+ + ICH,IQS,I3S,ISI,I3C
+
C==================================================================
C---fm to 1/GeV
DO 3 J1=1,30
RD(JJ)=RDH(LL,JJ)
EB(JJ)=EBH(LL,JJ)
RB(JJ)=RBH(LL,JJ)
+ IF(LL.NE.8.AND.LL.NE.9)GOTO25
C---Resets FD and RD for a nucleon-deuteron system (LL=8,9)
- IF(LL.EQ.8.OR.LL.EQ.9)THEN
- JH=LL-7+2*JJ-2
- FD(JJ)=AAND(1,JH)
- RD(JJ)=AAND(2,JH)-2*AAND(3,JH)/AAND(1,JH)
- ENDIF
+ JH=LL-7+2*JJ-2
+ FD(JJ)=AAND(1,JH)
+ RD(JJ)=AAND(2,JH)-2*AAND(3,JH)/AAND(1,JH)
C---Resets FD and RD for a pion-pion system (LL=5,6,7)
- IF(LL.EQ.5.OR.LL.EQ.6.OR.LL.EQ.7)THEN
- IF(LL.EQ.7)FD(JJ)=AAPI(1,2)/AM
- IF(LL.EQ.5)FD(JJ)=(.6667D0*AAPI(1,1)+.3333D0*AAPI(1,2))/AM
- IF(LL.EQ.6)FD(JJ)=(.3333D0*AAPI(1,1)+.6667D0*AAPI(1,2))/AM
- AKS=0.D0
- DAKS=1.D-5
- AKSH=AKS+DAKS
- AKH=DSQRT(AKSH)
- GPI1H=GPIPI(AKSH,1)
- GPI2H=GPIPI(AKSH,2)
- H=1/FD(JJ)
- IF(LL.EQ.7)C(JJ)=1/DCMPLX(GPI2H,-AKH)
- IF(LL.EQ.5)
- + C(JJ)=.6667D0/DCMPLX(GPI1H,-AKH)+.3333D0/DCMPLX(GPI2H,-AKH)
- IF(LL.EQ.6)
- + C(JJ)=.3333D0/DCMPLX(GPI1H,-AKH)+.6667D0/DCMPLX(GPI2H,-AKH)
- HH=DREAL(1/C(JJ))
- RD(JJ)=2*(HH-H)/DAKS
- ENDIF
-C---Resets FD and RD for a pion-nucleon system (LL=12,13)
- IF(LL.EQ.12.OR.LL.EQ.13)THEN
- IF(LL.EQ.12)FD(JJ)=AAPIN(1,2)
- IF(LL.EQ.13)FD(JJ)=(.6667D0*AAPIN(1,1)+.3333D0*AAPIN(1,2))
- AKS=0.D0
- DAKS=1.D-5
- AKSH=AKS+DAKS
- AKH=DSQRT(AKSH)
- GPI1H=GPIN(AKSH,1)
- GPI2H=GPIN(AKSH,2)
- H=1/FD(JJ)
- IF(LL.EQ.12)C(JJ)=1/DCMPLX(GPI2H,-AKH)
- IF(LL.EQ.13)
- + C(JJ)=.6667D0/DCMPLX(GPI1H,-AKH)+.3333D0/DCMPLX(GPI2H,-AKH)
- HH=DREAL(1/C(JJ))
- RD(JJ)=2*(HH-H)/DAKS
- ENDIF
+ 25 IF(LL.GT.7.OR.LL.LT.5)GOTO 24
+ IF(LL.EQ.7)FD(JJ)=AAPI(1,2)/AM
+ IF(LL.EQ.5)FD(JJ)=(.6667D0*AAPI(1,1)+.3333D0*AAPI(1,2))/AM
+ IF(LL.EQ.6)FD(JJ)=(.3333D0*AAPI(1,1)+.6667D0*AAPI(1,2))/AM
+ AKS=0.D0
+ DAKS=1.D-5
+ AKSH=AKS+DAKS
+ AKH=DSQRT(AKSH)
+ GPI1H=GPIPI(AKSH,1)
+ GPI2H=GPIPI(AKSH,2)
+ H=1/FD(JJ)
+ IF(LL.EQ.7)C(JJ)=1/DCMPLX(GPI2H,-AKH)
+ IF(LL.EQ.5)
+ +C(JJ)=.6667D0/DCMPLX(GPI1H,-AKH)+.3333D0/DCMPLX(GPI2H,-AKH)
+ IF(LL.EQ.6)
+ +C(JJ)=.3333D0/DCMPLX(GPI1H,-AKH)+.6667D0/DCMPLX(GPI2H,-AKH)
+ HH=DREAL(1/C(JJ))
+ RD(JJ)=2*(HH-H)/DAKS
+ 24 CONTINUE
C---Calculation continues for any system (any LL)
55 CONTINUE
- RETURN
END
+C=======================================================
+
+
+
+
+
+
+
+
+
- FUNCTION GPIPI(X,J)
+ FUNCTION GPIPIold(X,J)
C--- GPIPI = k*COTG(DELTA), X=k^2
C-- J=1(2) corresponds to isospin=0(2)
IMPLICIT REAL*8 (A-H,O-Z)
1 X2,Y2,Z2,T2,R2
COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
COMMON/FSI_FFPN/FF12,FF21
- COMPLEX*16 FF12,FF21
COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS
-
+ COMPLEX*16 FF12,FF21
C------------------------------------------------------------------
-c write(*,*)'in fsiw C1 C2 CN',C1,C2,CN
-
C==> AC1,2 = "relativistic" Bohr radii for particle-nucleus systems
-
C1N=C1*CN
IF(C1N.NE.0.D0)AC1=137.036D0/(C1N*E1) !m1-->E1
C2N=C2*CN
IF(C2N.NE.0.D0)AC2=137.036D0/(C2N*E2) !m2-->E2
C-----------------------------------------------------------
- CALL FSIPN !(WEIF) !weight due to particle-nucleus FSI
+ J=1
+ CALL FSIPN !weight due to particle-nucleus FSI
JRAT=0
- CALL FSIWF !(WEI) !weight due to particle-particle-nucleus FSI
+ CALL FSIWF !weight due to particle-particle-nucleus FSI
WEIN=WEI
- IF(I3C*J.NE.0) THEN
+ IF(I3C*J.NE.0) THEN
FF12=DCMPLX(1.D0,0.D0)
FF21=DCMPLX(1.D0,0.D0)
JRAT=1
- CALL VZ !(WEIN) ! weight due to particle-particle FSI
- ENDIF
+ CALL VZ ! weight due to particle-particle FSI
+ ENDIF
RETURN
END
+
C=======================================================================
SUBROUTINE LTRAN(P0,P,PS)
RETURN
END
- FUNCTION GND(X,J)
-C--- GND = k*COTG(DELTA), X=k^2
-C--- J=1(2) corresp. to nd(pd), S=1/2,
-C--- J=3(4) corresp. to nd(pd), S=3/2
- IMPLICIT REAL*8 (A-H,O-Z)
- COMMON/FSI_AAND/AAND(20,4)
- XX=X
- GND=1/AAND(1,J)+.5D0*AAND(2,J)*X
- DO 1 I=4,4
- XX=XX*X
- 1 GND=GND+AAND(I,J)*XX
- GND=GND/(1+AAND(3,J)*X)
- RETURN
- END
FUNCTION GDD(X,J)
C--- GDD = k*COTG(DELTA), X=k^2
RETURN
END
-
+
+ SUBROUTINE FSIPN
+C calculating particle-nucleus Coulomb Wave functions FFij
+ IMPLICIT REAL*8 (A-H,O-Z)
+ COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
+ COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1, !part. momenta in NRF
+ 1 P2X,P2Y,P2Z,E2,P2
+ COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, ! 4-coord. of emis. points in NRF
+ 1 X2,Y2,Z2,T2,R2
+ COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
+ COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
+ COMMON/FSI_ICH1/ICH1
+ COMMON/FSI_ETA/ETA
+ COMMON/FSI_FFPN/FF12,FF21
+ COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS
+
+ COMPLEX*16 FF1,FF12,FF21
+ FF12=DCMPLX(1.D0,0.D0)
+ FF21=DCMPLX(1.D0,0.D0)
+ IF(I3C.EQ.0)RETURN
+ ICH1=IDINT(C1)
+ IF(ICH1.EQ.0)GOTO 11
+ XH=AC1*P1
+ ACH=ACP(XH)
+ ACHR=DSQRT(ACH)
+ ETA=0.D0
+ IF(XH.NE.0.D0)ETA=1/XH
+ RHOS=P1*R1
+ HS=X1*P1X+Y1*P1Y+Z1*P1Z
+ FF12=FF12*FF1(RHOS,HS)
+ IF(IQS.EQ.0)GOTO 11
+ RHOS=P1*R2
+ HS=X2*P1X+Y2*P1Y+Z2*P1Z
+ FF21=FF21*FF1(RHOS,HS)
+ 11 ICH1=IDINT(C2)
+ IF(ICH1.EQ.0)GOTO 10
+ XH=AC2*P2
+ ACH=ACP(XH)
+ ACHR=DSQRT(ACH)
+ ETA=0.D0
+ IF(XH.NE.0.D0)ETA=1/XH
+ RHOS=P2*R2
+ HS=X2*P2X+Y2*P2Y+Z2*P2Z
+ FF12=FF12*FF1(RHOS,HS)
+CW WRITE(6,41)'AC2 ',AC2,'ACH ',ACH,'ETA ',ETA,'RHOS ',RHOS,'HS ',HS
+41 FORMAT(5(A5,E11.4))
+CW WRITE(6,40)'FF12 ',DREAL(FF12),DIMAG(FF12)
+ IF(IQS.EQ.0)GOTO 10
+ RHOS=P2*R1
+ HS=X1*P2X+Y1*P2Y+Z1*P2Z
+ FF21=FF21*FF1(RHOS,HS)
+CW WRITE(6,41)'AC1 ',AC1,'ACH ',ACH,'ETA ',ETA,'RHOS ',RHOS,'HS ',HS
+CW WRITE(6,40)'FF21 ',DREAL(FF21),DIMAG(FF21)
+40 FORMAT(A7,2E12.4)
+ 10 CONTINUE
+
+C WEIF = the weight due to the Coulomb particle-nucleus interaction
+ WEIF=DREAL(FF12)**2+DIMAG(FF12)**2
+ IF(IQS.EQ.1)WEIF=0.5D0*(WEIF+DREAL(FF21)**2+DIMAG(FF21)**2)
+ RETURN
+ END
+
+ FUNCTION GPIPI(X,J)
+C--- GPIPI = k*COTG(DELTA), X=k^2
+C-- J=1(2) corresponds to isospin=0(2)
+ IMPLICIT REAL*8 (A-H,O-Z)
+ COMMON/FSI_AAPI/AAPI(20,2)
+ COMMON/FSI_C/HELP(20),AM,AMS,DM
+ OM=DSQRT(X+AMS)
+ XX=X/AMS
+ GPIPI=OM/AAPI(1,J)
+ GPIPI=GPIPI*(1+(AAPI(3,J)-AAPI(1,J)**2)*XX+AAPI(4,J)*XX*XX)
+ GPIPI=GPIPI/(1+(AAPI(3,J)+AAPI(2,J)/AAPI(1,J))*XX)
+ RETURN
+ END
+ FUNCTION GND(X,J)
+C--- GND = k*COTG(DELTA), X=k^2
+C--- J=1(2) corresp. to nd(pd), S=1/2,
+C--- J=3(4) corresp. to nd(pd), S=3/2
+ IMPLICIT REAL*8 (A-H,O-Z)
+ COMMON/FSI_AAND/AAND(20,4)
+ XX=X
+ GND=1/AAND(1,J)+.5D0*AAND(2,J)*X
+ DO 1 I=4,4
+ XX=XX*X
+ 1 GND=GND+AAND(I,J)*XX
+ GND=GND/(1+AAND(3,J)*X)
+ RETURN
+ END
+
SUBROUTINE CKKB ! calculates KK-b scattering amplitude,
! saturated by S*(980) and delta(982) resonances
IMPLICIT REAL*8 (A-H,O-Z)
C(1)=C(1)+AAKK(8)/2/DCMPLX(AAKK(5)-S,
,-AK*AAKK(8)-AKPIETA*AAKK(9))
RETURN
- END
-
- SUBROUTINE VZ !(WEI)
+ END
+
+
+C
+ SUBROUTINE FSIWF
+C==> Prepares necessary quantities and call VZ(WEI) to calculate
+C the weight due to FSI
+ IMPLICIT REAL*8 (A-H,O-Z)
+ COMMON/FSI_CVK/V,CVK
+ COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1, !part. momenta in NRF
+ 1 P2X,P2Y,P2Z,E2,P2
+ COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS,
+ 1 X,Y,Z,T,RP,RPS
+ COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, ! 4-coord. of emis. points in NRF
+ 1 X2,Y2,Z2,T2,R2
+ COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
+ COMMON/FSI_SPIN/RHO(10)
+ COMMON/FSI_BP/B,P
+ COMMON/FSI_ETA/ETA
+ COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
+ COMMON/FSI_SW/RB(10),EB(10),BK(10),CDK(10),SDK(10),
+ 1 SBKRB(10),SDKK(10)
+ COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
+ COMMON/FSI_RR/F(10)
+ COMMON/FSI_FD/FD(10),RD(10)
+ COMMON/FSI_C/C(10),AM,AMS,DM
+ COMPLEX*16 C,F
+ COMMON/FSI_AA/AA
+ COMMON/FSI_SHH/SH,CHH
+ COMMON/FSI_AAPI/AAPI(20,2)/FSI_AAND/AAND(20,4)
+ COMMON/FSI_P12/P12X,P12Y,P12Z,E12,P12,AM12,EPM
+ COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS
+
+
+C==>calculating relative 4-coordinates of the particles in PRF
+C- {T,X,Y,Z} from the relative coordinates {TS,XS,YS,ZS} in NRF
+ XS=X1-X2
+ YS=Y1-Y2
+ ZS=Z1-Z2
+ TS=T1-T2
+ RS12=XS*P12X+YS*P12Y+ZS*P12Z
+ H1=(RS12/EPM-TS)/AM12
+
+c-mlv X=XS+P12X*H1
+c-mlv Y=YS+P12Y*H1
+c-mlv Z=ZS+P12Z*H1
+c-mlv T=(E12*TS-RS12)/AM12
+c-mlv RPS=X*X+Y*Y+Z*Z
+c-mlv RP=DSQRT(RPS)
+c WRITE(6,38)'RP ',RP,'X ',X,Y,Z,T
+c38 FORMAT(A7,E11.4,A7,4E11.4)
+
+ CVK=(P12X*PPX+P12Y*PPY+P12Z*PPZ)/(P12*AK)
+ V=P12/E12
+
+ IF(ICH.EQ.0)GOTO 21
+ XH=AC*AK
+ ACH=ACP(XH)
+ ACHR=DSQRT(ACH)
+ ETA=0.D0
+ IF(XH.NE.0.D0)ETA=1/XH
+C---HCP, HPR needed (e.g. in GST) if ICH=1
+ HCP=HC(XH)
+ HPR=HCP+.1544313298D0
+ 21 CONTINUE
+ MSP=MSPIN
+ DO 30 JJ=1,MSP
+ IF(NS.NE.1)GOTO22
+C---Calc. quantities for the square well potential;
+C-- for LL > 5 the square well potential is not possible or available
+cc GAK=G(AK)
+ BK(JJ)=DSQRT(EB(JJ)**2+AKS)
+ XRA=2*RB(JJ)/AC
+ HRA=BK(JJ)*RB(JJ)
+ CALL SEQ(XRA,HRA)
+ SBKRB(JJ)=HRA*B
+ HRA=AK*RB(JJ)
+ CALL GST(XRA,HRA)
+ SDK(JJ)=SH
+ CDK(JJ)=CHH
+ SDKK(JJ)=RB(JJ)
+ IF(AK.NE.0.D0)SDKK(JJ)=SH/AK
+ IF(ICH.EQ.1)SDK(JJ)=ACH*SDK(JJ)
+ 22 CONTINUE
+C-----------------------------------------------------------------------
+C---Calc. the strong s-wave scattering amplitude = C(JJ)
+C-- divided by Coulomb penetration factor squared (if ICH=1)
+ IF(NS.NE.1)GOTO 230
+ IF(LL.NE.4)GOTO 230 ! SW scat. amplitude used for alfa-alfa only
+ AKACH=AK
+ IF(ICH.EQ.1)AKACH=AK*ACH
+ C(JJ)=1/DCMPLX(GAK,-AKACH) ! amplitude for the SW-potential
+ GOTO 30
+ 230 IF(LL.EQ.5.OR.LL.EQ.6.OR.LL.EQ.7)GOTO20 ! pipi
+ IF(LL.EQ.8.OR.LL.EQ.9)GOTO20 ! Nd
+ IF(LL.EQ.14.OR.LL.EQ.17.OR.LL.EQ.23)GOTO27 ! K+K-, K-p, K0K0-b
+ A1=RD(JJ)*FD(JJ)*AKS
+ A2=1+.5D0*A1
+ IF(ICH.EQ.1)A2=A2-2*HCP*FD(JJ)/AC
+ AKF=AK*FD(JJ)
+ IF(ICH.EQ.1)AKF=AKF*ACH
+ C(JJ)=FD(JJ)/DCMPLX(A2,-AKF)
+ GOTO30
+ 20 CONTINUE
+C---Calc. scatt. ampl. C(JJ) for pipi and Nd
+ JH=LL-7+2*JJ-2
+ IF(LL.GT.7)GPI2=GND(AKS,JH)
+ IF(LL.LE.7)GPI2=GPIPI(AKS,2)
+ C(JJ)=1.D0/DCMPLX(GPI2,-AK)
+ IF(LL.GE.7)GOTO27
+ GPI1=GPIPI(AKS,1)
+ IF(LL.EQ.5)C(JJ)=.6667D0/DCMPLX(GPI1,-AK)+.3333D0*C(JJ)
+ IF(LL.EQ.6)C(JJ)=.3333D0/DCMPLX(GPI1,-AK)+.6667D0*C(JJ)
+ 27 CONTINUE
+C---Calc. K+K-, K0K0-b or K-p s-wave scatt. ampl.
+ IF(LL.EQ.14.OR.LL.EQ.23)CALL CKKB
+ IF(LL.EQ.17)C(JJ)=DCMPLX(3.29D0,3.55D0)
+C---Calc. pi+pi-, pi+pi+, pd, K+K- or K-p s-wave scatt. amplitude
+C-- divided by Coulomb penetration factor squared (if ICH=1)
+ IF(ICH.EQ.0)GOTO 30
+ AAK=ACH*AK
+ HCP2=2*HCP/AC
+ C(JJ)=1/(1/C(JJ)-HCP2+DCMPLX(0.D0,AK-AAK))
+ 30 CONTINUE
+C***********************************************************************
+ CALL VZ
+ RETURN
+ END
+
+ SUBROUTINE VZ
C==> Calculates the weight WEI due to FSI
IMPLICIT REAL*8 (A-H,O-Z)
COMMON/FSI_JR/JRAT
COMMON/FSI_C/C(10),AM,AMS,DM
COMMON/FSI_COULPH/EIDC
COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS
+
COMPLEX*16 F,C,G,PSI12,PSI21
COMPLEX*16 F12,F21
COMPLEX*16 EIDC
COMPLEX*8 Z8,CGAMMA
COMMON/FSI_FFPN/FF12,FF21
COMPLEX*16 FF12,FF21
-c--mlv------------------------------------------------------
-c write(*,*)'------- WE are in VZ-----------'
-c write(*,*)'JRAT=', JRAT
-c write(*,*)'AK=', AK,'RP=',RP
-c write(*,*)'X Y Z', X,Y,Z
-c write(*,*)'PPX PPY PPZ',PPX,PPY,PPZ
-c write(*,*)'ETA',ETA
-c write(*,*)'MSPIN',MSPIN
-c write(*,*)'ISI',ISI
-c write(*,*)'RP RPS',RP,RPS
-c write(*,*)'AA- BEFORE FIRT------',AA
-c write(*,*)'IQS',IQS
-c write(*,*)'ICH',ICH
-c write(*,*)'MSPIN',MSPIN
-c write(*,*)'C',(C(I),I=1,MSPIN)
-c write(*,*)'ACHR',ACHR
-c write(*,*)'F12',F12,'F21',F21,'FF12',FF12,'FF21',FF21
-
-c--mlv------------------------------------------------------
WEI=0.D0
IF(JRAT.EQ.1)GOTO 11
RHOS=AK*RP
HS=X*PPX+Y*PPY+Z*PPZ
-c write(*,*)'rhos',rhos,'hs',hs
IF(RHOS.LT.15.D0.AND.RHOS+DABS(HS).LT.20.D0)GOTO 2
C---Calc. EIDC=exp(i*Coul.Ph.);
C-- used in calc. of hypergeom. f-s in SEQA, FAS at k*R > 15, 20
-c write(*,*)'1111'
Z8=CMPLX(1.,SNGL(ETA))
- Z8=CGAMMA(Z8)
+ Z8=CGAMMA(Z8)
EIDC=Z8/CABS(Z8)
-c write(*,*)'1111 z8 eidc',z8,eidc
+
+c write(*,*)'Z8 EIDC',Z8,EIDC
2 CALL FF(RHOS,HS)
11 MSP=MSPIN
-
-c write(*,*)'before go to 4', wei
-
IF(ISI.EQ.0)GOTO 4 ! the strong interaction ON (OFF) if ISI=1(0)
IF(RP.LT.AA)GOTO 4
IF(JRAT.NE.1) CALL FIRT
PSI12=FF12*(F12+G)
PSI21=FF21*(F21+G)
G=PSI12+JSIGN*PSI21
- 1 WEI=WEI+RHO(JJ)*(DREAL(G)**2+DIMAG(G)**2)
-
-c WRITE(*,*)'WEI111=',WEI
-
+ 1 WEI=WEI+RHO(JJ)*(DREAL(G)**2+DIMAG(G)**2)
GOTO 8
5 DO 6 JJ=1,MSP
G=F(JJ)*C(JJ)
PSI12=FF12*(F12+G)
6 WEI=WEI+RHO(JJ)*(DREAL(PSI12)**2+DIMAG(PSI12)**2)
RETURN
-
-c write(*,*)'before 4 wei',wei
4 PSI12=FF12*F12
IF(IQS.EQ.0)GOTO 50 ! the quantum statistics ON (OFF) if IQS=1(0)
PSI21=FF21*F21
DO 3 JJ=1,MSP
JSIGN=-JSIGN
G=PSI12+JSIGN*PSI21
-c write(*,*)'msp jsign psi12 psi21',msp,jsign,psi12,psi21
3 WEI=WEI+RHO(JJ)*(DREAL(G)**2+DIMAG(G)**2)
-c write(*,*)'wei rho(1)G dreal(G) dimag(G)',wei,rho(1),G,
-c + dreal(G),dimag(G)
GOTO 8
-
50 WEI=DREAL(PSI12)**2+DIMAG(PSI12)**2
-c write(*,*)'WEI ---in VZ 50 before returne', WEI
RETURN
8 WEI=WEI/2
-c write(*,*)'WEI ---in VZ 8 before returne and END', WEI
+
+c [M %/WRITE(*,*)WEI
+
RETURN
END
-
+
+
SUBROUTINE FIRT
C---CALC. THE F(JJ)
C-- F(JJ)*C(JJ)= DEVIATION OF THE BETHE-SALPETER AMPL. FROM PLANE WAVE
IF(AK.NE.0.D0)GOTO4
HM=H*RPS
IF(HM.GE.3.D15)GOTO6
-
FS1=DFRSIN(HM)
FC1=DFRCOS(HM)
-
FC2=FC1
FS2=FS1
GOTO5
HM=H*(RSS-H1)**2
HP=H*(RSS+H1)**2
IF(HP.GE.3.D15)GOTO6
-
FS1=DFRSIN(HM)
FC1=DFRCOS(HM)
FS2=DFRSIN(HP)
FC2=DFRCOS(HP)
-
GOTO 5
11 CONTINUE
FS1=0.D0
IF(AK.NE.0.D0)GOTO 12
HM=H*RPS
IF(HM.GE.3.D15)GOTO6
-
FS1=FS1+DFRSIN(HM)/2
FC1=FC1+DFRCOS(HM)/2
-
IF(I.EQ.1)GOTO 13
FC2=FC1
FS2=FS1
HM=H*(RSS-H1)**2
HP=H*(RSS+H1)**2
IF(HP.GE.3.D15)GOTO6
-
FS1=FS1+DFRSIN(HM)/2
FC1=FC1+DFRCOS(HM)/2
FS2=FS2+DFRSIN(HP)/2
FC2=FC2+DFRCOS(HP)/2
-
13 CONTINUE
5 C12=FC1+FS2
S12=FS1+FC2
10 CONTINUE
RETURN
END
-
FUNCTION EXF(X)
IMPLICIT REAL*8 (A-H,O-Z)
IF(X.LT.-15.D0) GO TO 1
1 EXF=.0D0
RETURN
END
-
SUBROUTINE SEQ(X,H)
C---CALC. FUNCTIONS B, P (EQS. (17) OF G-K-L-L);
C-- NEEDED TO CALC. THE CONFLUENT HYPERGEOMETRIC FUNCTION GST.
PH(2)=PH(3)
GOTO 2
END
-
SUBROUTINE SEQA(X,H)
C---CALC. FUNCTIONS CHH=REAL(GST), SH=IMAG(GST)/ACH, B=SH/H
C-- IN THE ASYMPTOTIC REGION H=K*R >> 1.
B=SH/H
RETURN
END
-
SUBROUTINE FF(RHO,H)
C---Calc. F12, F21;
C-- F12= FF0* plane wave, FF0=F*ACHR,
F21=FF0(RHO,-H)*F21
RETURN
END
-
FUNCTION FAS(RKS)
C-- FAS=F*ACHR
C---F is the confluent hypergeometric function at k*r >> 1
FAS=FAS-DCMPLX(DCOS(RKS),DSIN(RKS))*ETA/RKS/ZZ1
RETURN
END
-
FUNCTION FF0(RHO,H)
C-- FF0=F*ACHR
C-- F is the confluent hypergeometric function
FF0=FF0*ACHR
RETURN
END
-
FUNCTION HC(XA)
C---HC = h-function of Landau-Lifshitz: h(x)=Re[psi(1-i/x)]+ln(x)
C-- psi(x) is the digamma function (the logarithmic derivative of
4 XP=X2*XP
RETURN
END
-
FUNCTION ACP(X)
C--- ACP = COULOMB PENETRATION FACTOR
IMPLICIT REAL*8 (A-H,O-Z)
1 ACP=1.D-6
RETURN
END
-
SUBROUTINE GST(X,H)
C---CALC. THE CONFL. HYPERGEOM. F-N = CHH+i*SH
C-- AND THE COULOMB F-S B, P (CALLS SEQ OR SEQA).
4 CALL SEQA(X,H)
RETURN
END
-
+
FUNCTION FF1(RHO,H)
C---FF1=FF0; used for particle-nucleus system
C-- FF0=F12*ACHR
IF(RHO.LT.15.D0.AND.RHO+H.LT.20.D0)GOTO 2
C---Calc. EIDC=exp(i*Coul.Ph.);
C-- used in calc. of hypergeom. f-s in SEQA, FAS at k*R > 15, 20
-
Z8=CMPLX(1.,SNGL(ETA))
Z8=CGAMMA(Z8)
-
EIDC=Z8/CABS(Z8)
-
2 FF1=FF0(RHO,H)
RETURN
END
FUNCTION G(AK)
C---Used to calculate SW scattering amplitude for alpa-alpha system
-C-- and for sear.f (square well potential search)
C---NOTE THAT SCATT. AMPL.= 1/CMPLX(G(AK),-AK*ACH)
IMPLICIT REAL*8 (A-H,O-Z)
COMMON/FSI_SW/RB(10),EB(10),BK(10),CDK(10),SDK(10),
COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
COMMON/FSI_ACH/HPR,AC,ACH,ACHR,JJ,MSPIN
COMMON/FSI_BP/B,P/FSI_DERIV/BPR,PPR/FSI_SHH/SH,CHH
- COMMON/FSI_DAK/DAK,HAC,IFUN
HAC=.0D0
IF(ICH.EQ.0)GOTO 1
XH=AC*AK
C G= GG/[B(kp,d)*B'(k,d)-B'(kp,d)*B(k,d)]
RETURN
END
-
SUBROUTINE DERIW(X,H)
C---CALLED BY F-N G(AK)
IMPLICIT REAL*8 (A-H,O-Z)
PPR=PPR+X*(P-Q1)/HHH
RETURN
END
+C================================================================
+
c-------------#include "gen/pilot.h"
FUNCTION DFRSIN(X)
- IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+
+ IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION A(0:16),B(0:15),C1(0:25),C2(0:28)
+
PARAMETER (Z1 = 1, R8 = Z1/8, R32 = Z1/32)
+
DATA C0 /1.25331 41373 15500 3D0/
+
DATA NA,NB,NC1,NC2 /16,15,25,28/
+
DATA A( 0) / 0.76435 13866 41860 002D0/
DATA A( 1) /-0.43135 54754 76601 793D0/
DATA A( 2) / 0.43288 19997 97266 531D0/
END
c--#include "gen/pilot.h"
- FUNCTION CGAMMA(Z)
+ FUNCTION CGAMMA(Z)
+
COMPLEX*8 CGAMMA
COMPLEX*8 Z,U,V,F,H,S
CHARACTER NAME*(*)
CHARACTER*80 ERRTXT
PARAMETER (NAME = 'CGAMMA')
+
DIMENSION C(0:15)
+
PARAMETER (Z1 = 1, HF = Z1/2)
c--#if defined(CERNLIB_QF2C)
101 FORMAT('ARGUMENT EQUALS NON-POSITIVE INTEGER = ',1P,E15.1)
END
-
- SUBROUTINE FSIPN !(WEIF)
-C calculating particle-nucleus Coulomb Wave functions FFij
- IMPLICIT REAL*8 (A-H,O-Z)
- COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
- COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1, !part. momenta in NRF
- 1 P2X,P2Y,P2Z,E2,P2
- COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, ! 4-coord. of emis. points in NRF
- 1 X2,Y2,Z2,T2,R2
- COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
- COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
- COMMON/FSI_ICH1/ICH1
- COMMON/FSI_ETA/ETA
- COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS
- COMMON/FSI_FFPN/FF12,FF21
- COMPLEX*16 FF1,FF12,FF21
- FF12=DCMPLX(1.D0,0.D0)
- FF21=DCMPLX(1.D0,0.D0)
- IF(I3C.EQ.0)RETURN
- ICH1=IDINT(C1)
- IF(ICH1.EQ.0)GOTO 11
- XH=AC1*P1
- ACH=ACP(XH)
- ACHR=DSQRT(ACH)
- ETA=0.D0
- IF(XH.NE.0.D0)ETA=1/XH
- RHOS=P1*R1
- HS=X1*P1X+Y1*P1Y+Z1*P1Z
- FF12=FF12*FF1(RHOS,HS)
- IF(IQS.EQ.0)GOTO 11
- RHOS=P1*R2
- HS=X2*P1X+Y2*P1Y+Z2*P1Z
- FF21=FF21*FF1(RHOS,HS)
- 11 ICH1=IDINT(C2)
- IF(ICH1.EQ.0)GOTO 10
- XH=AC2*P2
- ACH=ACP(XH)
- ACHR=DSQRT(ACH)
- ETA=0.D0
- IF(XH.NE.0.D0)ETA=1/XH
- RHOS=P2*R2
- HS=X2*P2X+Y2*P2Y+Z2*P2Z
- FF12=FF12*FF1(RHOS,HS)
-CW WRITE(6,41)'AC2 ',AC2,'ACH ',ACH,'ETA ',ETA,'RHOS ',RHOS,'HS ',HS
-41 FORMAT(5(A5,E11.4))
-CW WRITE(6,40)'FF12 ',DREAL(FF12),DIMAG(FF12)
- IF(IQS.EQ.0)GOTO 10
- RHOS=P2*R1
- HS=X1*P2X+Y1*P2Y+Z1*P2Z
- FF21=FF21*FF1(RHOS,HS)
-CW WRITE(6,41)'AC1 ',AC1,'ACH ',ACH,'ETA ',ETA,'RHOS ',RHOS,'HS ',HS
-CW WRITE(6,40)'FF21 ',DREAL(FF21),DIMAG(FF21)
-40 FORMAT(A7,2E12.4)
- 10 CONTINUE
-
-C WEIF = the weight due to the Coulomb particle-nucleus interaction
- WEIF=DREAL(FF12)**2+DIMAG(FF12)**2
- IF(IQS.EQ.1)WEIF=0.5D0*(WEIF+DREAL(FF21)**2+DIMAG(FF21)**2)
- RETURN
- END
-
-C=======================================================
-C
- SUBROUTINE FSIWF !(WEI)
-C==> Prepares necessary quantities and call VZ(WEI) to calculate
-C the weight due to FSI
- IMPLICIT REAL*8 (A-H,O-Z)
- COMMON/FSI_CVK/V,CVK
- COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1, !part. momenta in NRF
- 1 P2X,P2Y,P2Z,E2,P2
- COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS,
- 1 X,Y,Z,T,RP,RPS
- COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, ! 4-coord. of emis. points in NRF
- 1 X2,Y2,Z2,T2,R2
- COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
- COMMON/FSI_SPIN/RHO(10)
- COMMON/FSI_BP/B,P
- COMMON/FSI_ETA/ETA
- COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
- COMMON/FSI_SW/RB(10),EB(10),BK(10),CDK(10),SDK(10),
- 1 SBKRB(10),SDKK(10)
- COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
- COMMON/FSI_RR/F(10)
- COMMON/FSI_FD/FD(10),RD(10)
- COMMON/FSI_C/C(10),AM,AMS,DM
- COMPLEX*16 C,F
- COMMON/FSI_AA/AA
- COMMON/FSI_SHH/SH,CHH
- COMMON/FSI_AAPI/AAPI(20,2)/FSI_AAND/AAND(20,4)
- COMMON/FSI_AAPIN/AAPIN(20,2)
- COMMON/FSI_P12/P12X,P12Y,P12Z,E12,P12,AM12,EPM
- COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS
-
-
-cmlv
- IF(IRANPOS.EQ.0)THEN
-C==>calculating relative 4-coordinates of the particles in PRF
-C- {T,X,Y,Z} from the relative coordinates {TS,XS,YS,ZS} in NRF
- XS=X1-X2
- YS=Y1-Y2
- ZS=Z1-Z2
- TS=T1-T2
- RS12=XS*P12X+YS*P12Y+ZS*P12Z
- H1=(RS12/EPM-TS)/AM12
- X=XS+P12X*H1
- Y=YS+P12Y*H1
- Z=ZS+P12Z*H1
- T=(E12*TS-RS12)/AM12
- RPS=X*X+Y*Y+Z*Z
- RP=DSQRT(RPS)
-c WRITE(6,38)'RP ',RP,'X ',X,Y,Z,T
-c FORMAT(A7,E11.4,A7,4E11.4)
- ENDIF
-
- CVK=(P12X*PPX+P12Y*PPY+P12Z*PPZ)/(P12*AK)
- V=P12/E12
-
- IF(ICH.EQ.0)GOTO 21
- XH=AC*AK
- ACH=ACP(XH)
- ACHR=DSQRT(ACH)
- ETA=0.D0
- IF(XH.NE.0.D0)ETA=1/XH
-C---HCP, HPR needed (e.g. in GST) if ICH=1
- HCP=HC(XH)
- HPR=HCP+.1544313298D0
- 21 CONTINUE
- MSP=MSPIN
- DO 30 JJ=1,MSP
- ISPIN=JJ
- IF(NS.NE.1)GOTO22
-C---Calc. quantities for the square well potential;
-C-- for LL > 5 the square well potential is not possible or available
- IF(LL.EQ.4)GOTO 22
- BK(JJ)=DSQRT(EB(JJ)**2+AKS)
- XRA=2*RB(JJ)/AC
- HRA=BK(JJ)*RB(JJ)
- CALL SEQ(XRA,HRA)
- SBKRB(JJ)=HRA*B
- HRA=AK*RB(JJ)
- CALL GST(XRA,HRA)
- SDK(JJ)=SH
- CDK(JJ)=CHH
- SDKK(JJ)=RB(JJ)
- IF(AK.NE.0.D0)SDKK(JJ)=SH/AK
- IF(ICH.EQ.1)SDK(JJ)=ACH*SDK(JJ)
- 22 CONTINUE
-C-----------------------------------------------------------------------
-C---Calc. the strong s-wave scattering amplitude = C(JJ)
-C-- divided by Coulomb penetration factor squared (if ICH=1)
- IF(NS.NE.1)GOTO 230
- IF(LL.NE.4)GOTO 230 ! SW scat. amplitude used for alfa-alfa only
- GAK=G(AK)
- AKACH=AK
- IF(ICH.EQ.1)AKACH=AK*ACH
- C(JJ)=1/DCMPLX(GAK,-AKACH) ! amplitude for the SW-potential
- GOTO 30
- 230 IF(LL.EQ.5.OR.LL.EQ.6.OR.LL.EQ.7)GOTO20 ! pipi
- IF(LL.EQ.12.OR.LL.EQ.13)GOTO20 ! piN
- IF(LL.EQ.8.OR.LL.EQ.9.OR.LL.EQ.18)GOTO20 ! Nd, dd
- IF(LL.EQ.14.OR.LL.EQ.17.OR.LL.EQ.23)GOTO27 ! K+K-, K-p, K0K0-b
- A1=RD(JJ)*FD(JJ)*AKS
- A2=1+.5D0*A1
- IF(ICH.EQ.1)A2=A2-2*HCP*FD(JJ)/AC
- AKF=AK*FD(JJ)
- IF(ICH.EQ.1)AKF=AKF*ACH
- C(JJ)=FD(JJ)/DCMPLX(A2,-AKF)
- GOTO30
- 20 CONTINUE
-C---Calc. scatt. ampl. C(JJ) for pipi, piN and Nd, dd
- JH=LL-7+2*JJ-2
- IF(LL.EQ.8.OR.LL.EQ.9)GPI2=GND(AKS,JH)
- IF(LL.EQ.18)GPI2=GDD(AKS,JJ)
- IF(LL.EQ.5.OR.LL.EQ.6.OR.LL.EQ.7)GPI2=GPIPI(AKS,2)
- IF(LL.EQ.12.OR.LL.EQ.13)GPI2=GPIN(AKS,2)
- C(JJ)=1.D0/DCMPLX(GPI2,-AK) !pi+pi+, nd, pd, pi+p, dd
- IF(LL.NE.5.AND.LL.NE.6.AND.LL.NE.13)GOTO27
- IF(LL.EQ.5.OR.LL.EQ.6)GPI1=GPIPI(AKS,1)
- IF(LL.EQ.13)GPI1=GPIN(AKS,1)
- IF(LL.EQ.5.OR.LL.EQ.13)
- c C(JJ)=.6667D0/DCMPLX(GPI1,-AK)+.3333D0*C(JJ) !pi+pi-,pi-p
- IF(LL.EQ.6)C(JJ)=.3333D0/DCMPLX(GPI1,-AK)+.6667D0*C(JJ) !pi0pi0
- 27 CONTINUE
-C---Calc. K+K-, K0K0-b or K-p s-wave scatt. ampl.
- IF(LL.EQ.14.OR.LL.EQ.23)CALL CKKB
- IF(LL.EQ.17)C(JJ)=DCMPLX(3.29D0,3.55D0)
-C---Calc. pi+pi-, pi+pi+, pd, pi+p, pi-p, K+K- or K-p s-wave scatt. ampl.
-C-- divided by Coulomb penetration factor squared (if ICH=1)
- IF(ICH.EQ.0)GOTO 30
- AAK=ACH*AK
- HCP2=2*HCP/AC
- C(JJ)=1/(1/C(JJ)-HCP2+DCMPLX(0.D0,AK-AAK))
-c write(*,*)'c(jj)',c(jj)
- 30 CONTINUE
-C***********************************************************************
-c write(*,*)'before call vz in fsiwf wei',wei
- CALL VZ !(WEI)
- RETURN
- END
-
- SUBROUTINE ltran12
-C==>calculating particle momentum in PRF {EE,PPX,PPY,PPZ} from
-C- the momentum of the first particle {E1,P1X,P1Y,P1Z) in NRF
+
+ SUBROUTINE ltran12
+C
+
+c SUBROUTINE TRANS(icrf,irot)
+C==> TRANSformation to the Co-moving frame (icrf>0) and
+C Rotation to the system where (Pt || X),(irot=1).
+C-FSI ***************************************************
IMPLICIT REAL*8 (A-H,O-Z)
- COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1, !part. momenta in NRF
+ COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1, ! momenta in NRF
1 P2X,P2Y,P2Z,E2,P2
- COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS,
+ COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, ! 4-coord. of emis.
+ 1 X2,Y2,Z2,T2,R2 ! points in NRF
+ COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS, ! momenta in PRF
1 X,Y,Z,T,RP,RPS
COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
COMMON/FSI_P12/P12X,P12Y,P12Z,E12,P12,AM12,EPM
- COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, !4-coord. of emis. points in NRF
- 1 X2,Y2,Z2,T2,R2
- COMMON/FSI_CONS/PI,PI2,SPI,DR,W
- COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS
+ COMMON/FSI_CVK/V,CVK
+C-FSI ***************************************************
+ COMMON /PAIR/P12T,V12Z,GAMZ,V12T,CPHI,SPHI
-
-C fm --> 1/GeV
-c write(*,*)'in LTRAN12 W', W
-c write(*,*)'in LTRAN12 p1x p1y p1z', P1X,P1Y,P1Z
-c write(*,*)'in LTRAN12 p2x p2y p2z', P2X,P2Y,P2Z
-c write(*,*)'in LTRAN12 x y z t', X,Y,Z,T
-c write(*,*)'in LTRAN12 x2 y2 z2', X2,Y2,Z2,T2
-
-CMLV
- IF(IRANPOS.EQ.0)THEN
- X1=X1*W
- Y1=Y1*W
- Z1=Z1*W
- T1=T1*W
- X2=X2*W
- Y2=Y2*W
- Z2=Z2*W
- T2=T2*W
-C calculating Ri, Pi and Ei
- R1=DSQRT(X1*X1+Y1*Y1+Z1*Z1)
- R2=DSQRT(X2*X2+Y2*Y2+Z2*Z2)
- ENDIF
-CMLV
-
+ icrf=1
+ irot=1
+
+C---> Particle energies ---------
P1S=P1X*P1X+P1Y*P1Y+P1Z*P1Z
P2S=P2X*P2X+P2Y*P2Y+P2Z*P2Z
- P1=DSQRT(P1S)
- P2=DSQRT(P2S)
E1=DSQRT(AM1*AM1+P1S)
E2=DSQRT(AM2*AM2+P2S)
-C-----------------------------------------------------------------------
+C---> Pair parameters -----------
+ E12=E1+E2 ! Energy
+ P12X=P1X+P2X ! Px
+ P12Y=P1Y+P2Y ! Py
+ P12Z=P1Z+P2Z ! Pz
+ P12S=P12X**2+P12Y**2+P12Z**2
+ P12 =DSQRT(P12S)! Momentum
+ V12 =P12/E12 ! Velocity
+ CTH =P12Z/P12 ! cos(theta)
+ STH =DSQRT(1.D0-CTH**2) !sin
+ V12Z=V12*CTH ! Longit. V
+ GAMZ=1.D0/DSQRT(1.D0-V12Z**2)
+C-- V12T=V12*STH ! Transv. V in CMS (not needed)
+ P12TS=P12X*P12X+P12Y*P12Y
+ P12T=DSQRT(P12TS) !Pt
+C===> Azimuthal rotation (Pt||X) ============
+ IF(V12T.NE.0.D0) THEN
+ CPHI=P12X/P12T ! cos(phi)
+ SPHI=P12Y/P12T ! sin(phi)
+ IF((irot.eq.1)) THEN
+ CALL ROT8(P1X,P1Y,SPHI,CPHI,P1X,P1Y)
+ CALL ROT8(P2X,P2Y,SPHI,CPHI,P2X,P2Y)
+ CALL ROT8(X1,Y1,SPHI,CPHI,X1,Y1)
+ CALL ROT8(X2,Y2,SPHI,CPHI,X2,Y2)
+ END IF
+ ELSE ! Rotation impossible
+ CPHI=2.D0 ! to avoid
+ SPHI=2.D0 ! using it !
+ END IF
+C===> Co-moving ref. frame ============
+ IF(icrf.gt.0) THEN
+ CALL LTR8(P1Z,E1,V12Z,GAMZ,P1Z,E1a)
+ CALL LTR8(P2Z,E2,V12Z,GAMZ,P2Z,E2a)
+ P1S=P1X*P1X+P1Y*P1Y+P1Z*P1Z
+ P2S=P2X*P2X+P2Y*P2Y+P2Z*P2Z
+ E1=DSQRT(AM1*AM1+P1S)
+ E2=DSQRT(AM2*AM2+P2S)
+ CALL LTR8(Z1,T1,V12Z,GAMZ,Z1,T1)
+ CALL LTR8(Z2,T2,V12Z,GAMZ,Z2,T2)
+ END IF
+C===> Pair reference frame ============
+ P1=DSQRT(P1S)
+ P2=DSQRT(P2S)
E12=E1+E2
P12X=P1X+P2X
P12Y=P1Y+P2Y
P12Z=P1Z+P2Z
P12S=P12X**2+P12Y**2+P12Z**2
- AM12=DSQRT(E12**2-P12S)
+ P12 =DSQRT(P12S)
+ AM12S=E12*E12-P12S
+ AM12=DSQRT(AM12S)
EPM=E12+AM12
- P12=DSQRT(P12S)
P112=P1X*P12X+P1Y*P12Y+P1Z*P12Z
H1=(P112/EPM-E1)/AM12
PPX=P1X+P12X*H1
EE=(E12*E1-P112)/AM12
AKS=EE**2-AM1**2
AK=DSQRT(AKS)
-CW WRITE(6,38)'AK ',AK,'K ',PPX,PPY,PPZ,EE
-38 FORMAT(A7,E11.4,A7,4E11.4)
+ CVK=(P12X*PPX+P12Y*PPY+P12Z*PPZ)/(P12*AK)
+ V=P12/E12
+ V12T=P12T/SQRT(AM12S+P12TS) ! transverse velocity in LCMS
+C---> Coordinates -----------------------------
+ XS=X1-X2
+ YS=Y1-Y2
+ ZS=Z1-Z2
+ TS=T1-T2
+ RS12=XS*P12X+YS*P12Y+ZS*P12Z
+ H1=(RS12/EPM-TS)/AM12
+
+cmlv X=XS+P12X*H1
+cmlv Y=YS+P12Y*H1
+cmlv Z=ZS+P12Z*H1
+cmlv T=(E12*TS-RS12)/AM12
+cmlv RPS=X*X+Y*Y+Z*Z
+cmlv RP=DSQRT(RPS)
+
+
+ RETURN
+ END
+C====
+ SUBROUTINE LTR8(Z,T,BETA,GAMMA,ZT,TT)
+C===> Lorentz Transf. of Z(Pz) and T(E) to moving ref. frame.(REAL*8)
+CInp: Z,T-Zcoord,Time before tr., BETA,GAMMA- velocity, Lor.fact.
+COut: ZT,TT- " " after transformation.
+C==== ===============================================================
+ IMPLICIT REAL*8 (A-H,O-Z)
+ ZH=GAMMA*(Z-BETA*T)
+ TT=GAMMA*(T-BETA*Z)
+ ZT=ZH
RETURN
END
+C====
+ SUBROUTINE LTR4(Z,T,BETA,GAMMA,ZT,TT)
+C===> Lorentz Transf. of Z(Pz) and T(E) to moving ref. frame.(real*4)
+CInp: Z,T-Zcoord,Time before tr., BETA,GAMMA- velocity, Lor.fact.
+COut: ZT,TT- " " after transformation.
+C==== ===============================================================
+ ZH=GAMMA*(Z-BETA*T)
+ TT=GAMMA*(T-BETA*Z)
+ ZT=ZH
+ RETURN
+ END
+C====
+ SUBROUTINE ROT8(X,Y,SF,CF,XR,YR)
+C===> Rotation with the angle f. (REAL*8)
+CInp: X,Y-coord. before rotation; SF=sin(f), CF=cos(f),
+COut: XR,YR - coordinates after rotation.
+C==== =================================================
+ IMPLICIT REAL*8 (A-H,O-Z)
+ XH=X*CF+Y*SF !Y
+ YR=Y*CF-X*SF ! _-X'
+ XR=XH ! _- f
+ RETURN !------>
+ END ! X
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file
git://git.uio.no / u/mrichter/AliRoot.git / commitdiff