--- /dev/null
+c 1 2 3 4 5 6 7 8
+c---------|---------|---------|---------|---------|---------|---------|---------|
+*-- Author : R.Lednicky 20/01/95
+ SUBROUTINE FSIW(J,WEIF,WEI,WEIN)
+
+C=======================================================================
+C Calculates final state interaction (FSI) weights
+C WEIF = weight due to particle - (effective) nucleus FSI (p-N)
+C WEI = weight due to p-p-N FSI
+C WEIN = weight due to p-p FSI; note that WEIN=WEI if I3C=0;
+C note that if I3C=1 the calculation of
+C WEIN can be skipped by putting J=0
+C.......................................................................
+C Correlation Functions:
+C CF(p-p-N) = sum(WEI)/sum(WEIF)
+C CF(p-p) = sum(WEIN)/sum(1); here the nucleus is completely
+C inactive
+C CF(p-p-"N") = sum(WEIN*WEIF')/sum(WEIF'), where WEIN and WEIF'
+C are not correlated (calculated at different emission
+C points, e.g., for different events);
+C thus here the nucleus affects one-particle
+C spectra but not the correlation
+C.......................................................................
+C User must supply data file <fn> on unit NUNIT (e.g. =11) specifying
+C LL : particle pair
+C NS : approximation used to calculate Bethe-Salpeter amplitude
+C ITEST: test switch
+C If ITEST=1 then also following parameters are required
+C ICH : 1(0) Coulomb interaction between the two particles ON (OFF)
+C IQS : 1(0) quantum statistics for the two particles ON (OFF)
+C ISI : 1(0) strong interaction between the two particles ON (OFF)
+C I3C : 1(0) Coulomb interaction with residual nucleus ON (OFF)
+C This data file can contain other information useful for the user.
+C It is read by subroutines READINT4 and READREA8(4) (or READ_FILE).
+C----------------------------------------------------------------------
+C- LL 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
+C- part. 1: n p n a pi+ pi0 pi+ n p pi+ pi+ pi+ pi- K+ K+ K+ K-
+C- part. 2: n p p a pi- pi0 pi+ d d K- K+ p p K- K+ p p
+C NS=1 y/n: + + + + + - - - - - - - - - - - -
+C----------------------------------------------------------------------
+C- LL 18 19 20 21 22 23 24 25 26 27 28 29 30
+C- part. 1: d d t t K0 K0 d p p p n lam p
+C- part. 2: d a t a K0 K0b t t a lam lam lam pb
+C NS=1 y/n: - - - - - - - - - + + + -
+C----------------------------------------------------------------------
+C NS=1 Square well potential,
+C NS=3 not used
+C NS=4 scattered wave approximated by the spherical wave,
+C NS=2 same as NS=4 but the approx. of equal emission times in PRF
+C not required (t=0 approx. used in all other cases).
+C Note: if NS=2,4, the B-S amplitude diverges at zero distance r* in
+C the two-particle c.m.s.; user can specify a cutoff AA in
+C SUBROUTINE FSIINI, for example:
+C IF(NS.EQ.2.OR.NS.EQ.4)AA=5.D0 !! in 1/GeV --> AA=1. fm
+C---------------------------------------------------------------------
+C ITEST=1 any values of parameters ICH, IQS, ISI, I3C are allowed
+C and should be given in data file <fn>
+C ITEST=0 physical values of these parameters are put automatically
+C in FSIINI (their values are not required in data file)
+C=====================================================================
+C At the beginning of calculation user should call FSIINI,
+C which reads LL, NS, ITEST (and eventually ICH, IQS, ISI, I3C)
+C and initializes various parameters.
+C In particular the constants in
+C COMMON/FSI_CONS/PI,PI2,SPI,DR,W
+C may be useful for the user:
+C W=1/.1973D0 ! from fm to 1/GeV
+C PI=4*DATAN(1.D0)
+C PI2=2*PI
+C SPI=DSQRT(PI)
+C DR=180.D0/PI ! from radian to degree
+C _______________________________________________________
+C !! |Important note: all real quantities are assumed REAL*8 | !!
+C -------------------------------------------------------
+C For each event user should fill in the following information
+C in COMMONs (all COMMONs in FSI calculation start with FSI_):
+C ...................................................................
+C COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
+C Only
+C AMN = mass of the effective nucleus [GeV/c**2]
+C CN = charge of the effective nucleus [elem. charge units]
+C are required
+C ...................................................................
+C COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1, !part. momenta in the rest frame
+C 1 P2X,P2Y,P2Z,E2,P2 !of effective nucleus (NRF)
+C Only the components
+C PiX,PiY,PiZ [GeV/c]
+C in NRF are required.
+C To make the corresponding Lorentz transformation user can use the
+C subroutines LTRAN and LTRANB
+C ...................................................................
+C COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, ! 4-coord. of emission
+C 1 X2,Y2,Z2,T2,R2 ! points in NRF
+C The componets
+C Xi,Yi,Zi [fm]
+C and emission times
+C Ti [fm/c]
+C should be given in NRF with the origin assumed at the center
+C of the effective nucleus. If the effect of residual nucleus is
+C not calculated within FSIW, the NRF can be any fixed frame.
+C-----------------------------------------------------------------------
+C Before calling FSIW the user must call
+C CALL LTRAN12
+C Besides Lorentz transformation to pair rest frame:
+C (p1-p2)/2 --> k* it also transforms 4-coordinates of
+C emission points from fm to 1/GeV and calculates Ei,Pi and Ri.
+C Note that |k*|=AK in COMMON/FSI_PRF/
+C-----------------------------------------------------------------------
+C After making some additional filtering using k* (say k* < k*max)
+C or direction of vector k*,
+C user can finally call FSIW to calculate the FSI weights
+C to be used to construct the correlation function
+C=======================================================================
+
+ IMPLICIT REAL*8 (A-H,O-Z)
+ COMMON/FSI_JR/JRAT
+ COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
+ COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1, ! particle momenta in NRF
+ 1 P2X,P2Y,P2Z,E2,P2
+ COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS, ! k*=(p1-p2)/2 and x1-x2
+ 1 X,Y,Z,T,RP,RPS ! in pair rest frame (PRF)
+ 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_FFPN/FF12,FF21
+ COMPLEX*16 FF12,FF21
+C------------------------------------------------------------------
+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
+ JRAT=0
+ call BoostToPrf()
+ CALL FSIWF(WEI) !weight due to particle-particle-nucleus FSI
+ WEIN=WEI
+ 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
+ RETURN
+ END
+C=======================================================================
+
+ SUBROUTINE LTRAN(P0,P,PS)
+C==>calculating particle 4-momentum PS={PSX,PSY,PSZ,ES}
+C in rest frame of a system 0 with 4-momentum P0={P0X,P0Y,P0Z,E0}
+C from its 4-momentum P={PX,PY,PZ,E}
+
+ IMPLICIT REAL*8 (A-H,O-Z)
+ DIMENSION P0(4),P(4),PS(4)
+C-----------------------------------------------------------------------
+ P0S=P0(1)**2+P0(2)**2+P0(3)**2
+ AM0=DSQRT(P0(4)**2-P0S)
+ EPM=P0(4)+AM0
+ PP0=P(1)*P0(1)+P(2)*P0(2)+P(3)*P0(3)
+ H=(PP0/EPM-P(4))/AM0
+ PS(1)=P(1)+P0(1)*H
+ PS(2)=P(2)+P0(2)*H
+ PS(3)=P(3)+P0(3)*H
+ PS(4)=(P0(4)*P(4)-PP0)/AM0
+ RETURN
+ END
+
+ SUBROUTINE LTRANB(P0,PS,P)
+C==>calculating particle 4-momentum P={PX,PY,PZ,E}
+C from its 4-momentum PS={PSX,PSY,PSZ,ES}
+C in rest frame of a system 0 with 4-momentum P0={P0X,P0Y,P0Z,E0}
+
+ IMPLICIT REAL*8 (A-H,O-Z)
+ DIMENSION P0(4),P(4),PS(4)
+C-----------------------------------------------------------------------
+ P0S=P0(1)**2+P0(2)**2+P0(3)**2
+ AM0=DSQRT(P0(4)**2-P0S)
+ EPM=P0(4)+AM0
+ PSP0=PS(1)*P0(1)+PS(2)*P0(2)+PS(3)*P0(3)
+ HS=(PSP0/EPM+PS(4))/AM0
+ P(1)=PS(1)+P0(1)*HS
+ P(2)=PS(2)+P0(2)*HS
+ P(3)=PS(3)+P0(3)*HS
+ P(4)=(P0(4)*PS(4)+PSP0)/AM0
+ 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
+ IMPLICIT REAL*8 (A-H,O-Z)
+ 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_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
+
+C fm --> 1/GeV
+ 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)
+ 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-----------------------------------------------------------------------
+ E12=E1+E2
+ P12X=P1X+P2X
+ P12Y=P1Y+P2Y
+ P12Z=P1Z+P2Z
+ P12S=P12X**2+P12Y**2+P12Z**2
+ AM12=DSQRT(E12**2-P12S)
+ EPM=E12+AM12
+ P12=DSQRT(P12S)
+ P112=P1X*P12X+P1Y*P12Y+P1Z*P12Z
+ H1=(P112/EPM-E1)/AM12
+ PPX=P1X+P12X*H1
+ PPY=P1Y+P12Y*H1
+ PPZ=P1Z+P12Z*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)
+ RETURN
+ 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/FSI_FFPN/FF12,FF21
+ COMPLEX*16 FF1,FF12,FF21
+ FF12=DCMPLX(1.D0,0.D0)
+ FF21=DCMPLX(1.D0,0.D0)
+C ACH=1.D0
+C WEIF=1.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 GPIN(X,J)
+C--- GPIN = k*COTG(DELTA), X=k^2
+C-- J=1(2) corresponds to piN isospin=1/2(3/2)
+ IMPLICIT REAL*8 (A-H,O-Z)
+ COMMON/FSI_AAPIN/AAPIN(20,2)
+ GPIN=1/AAPIN(1,J)+.5D0*AAPIN(2,J)*X
+ 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
+C--- J=1,2,3 corresp. to S=0,1,2
+ IMPLICIT REAL*8 (A-H,O-Z)
+ COMMON/FSI_AADD/AADD(20,3)
+ COMMON/FSI_C/C(10),AM,AMS,DM
+ COMMON/FSI_CONS/PI,PI2,SPI,DR,W
+ COMPLEX*16 C
+ E=X/2/AM
+ ER=DSQRT(E)
+ IF(J.EQ.1)THEN
+ GDD=ER*(AADD(1,1)*DEXP(-E/AADD(2,1))-AADD(3,1))
+ GDD=GDD/DR ! from degree to radian
+ TAND=DTAN(GDD)
+ IF(TAND.EQ.0.D0)TAND=1.D-10
+ GDD=DSQRT(X)/TAND
+ END IF
+ IF(J.EQ.2)THEN
+ GDD=1.D10
+ END IF
+ IF(J.EQ.3)THEN
+ GDD=ER*(AADD(1,3)+AADD(2,3)*E)
+ GDD=GDD/DR ! from degree to radian
+ TAND=DTAN(GDD)
+ IF(TAND.EQ.0.D0)TAND=1.D-10
+ GDD=DSQRT(X)/TAND
+ END IF
+ RETURN
+ END
+ BLOCK DATA
+ IMPLICIT REAL*8 (A-H,O-Z)
+ COMMON/FSI_AAPI/AAPI(20,2)/FSI_AAND/AAND(20,4)
+ COMMON/FSI_AADD/AADD(20,3)/FSI_AAPIN/AAPIN(20,2)
+ COMMON/FSI_AAKK/AAKK(9)/FSI_AAPAP/AAPAPR(3,2),AAPAPI(3,2)
+C---Parameters for GPIPI (I,J), J=1,2 -> isospin=0,2
+ DATA AAPI/.2600D00, .2500D00, .3480D00,-.0637D00, 16*0.D0,
+ 1 -.0280D00,-.0820D00, .2795D00,-.0086D00, 16*0.D0/
+C---Parameters for GPIN (I,J), J=1,2 -> piN isospin=1/2,3/2
+ DATA AAPIN/ .12265D1, .1563D2, 18*0.D0,
+ 1 -.750D0, .7688D2, 18*0.D0/
+C---Parameters for GND (I,J), J=1(2) corresp. to nd(pd), S=1/2,
+C--- J=3(4) corresp. to nd(pd), S=3/2
+ DATA AAND/-.3295D1,-.8819D3, .4537D4, .28645D5, 16*0.D0,
+ 1 -.13837D2, .11505D2, .0D0, .10416D2, 16*0.D0,
+ 2 -.32180D2, .1014D2, .0D0, .0D0, 16*0.D0,
+ 3 -.60213D2, .1333D2, .0D0,-.70309D2, 16*0.D0/
+ DATA AADD/ .10617D4, .3194D-2, .56849D3, 17*0.D0,
+ 1 20*0.D0,
+ 2 -.1085D4, .1987D5, 18*0.D0/
+C--- AAKK= m_K^2, m_pi^2, m_eta^2, m_S*^2, m_delta^2,
+C--- gam(S*-->KK-b), gam(S*-->pipi), gam(delta-->KK-b),
+C--- gam(delta-->pi eta)
+ DATA AAKK/.247677D00,.01947977D00,.2997015D00,.9604D00,
+ 1 .96511D00,
+cc 2 .792D00, .199D00, .333D00, .222D00/ ! Martin (77)
+ 2 .094D00, .110D00, .333D00, .222D00/ ! Morgan (93)
+C---Parameters for PAP (I,J), j=1,2 -> isospin I=0,2
+C--- i=1-3 -> a_singlet, a_triplet, d [fm]
+C--- Im a_IS (I=isospin, S=spin) are fixed by atomic data and
+C n-bar survival time up to one free parameter, e.g. Im a_00
+C--- Batty (89), Kerbikov (93):
+C--- Ima_10=1.96-Ima_00, Ima_01=0.367-Ima_00/3, Ima_11=0.453+Ima_00/3
+C--- In DATA we put Ima_00=0.3.
+C--- Re a_IS are fixed by atomic data up to three free parameters
+C--- Batty (89):
+C--- Rea_aver(pp-bar)=Re[(a_00+a_10)+3(a_01+a_11)]/8=-0.9
+C--- In DATA we used Rea_IS from Paris potential Pignone (94)
+C--- rescaled by 1.67 to satisfy the atomic constraint.
+C--- Effective radius is is taken independent of IS from the phase
+C--- shift fit by Pirner et al. (91).
+ DATA AAPAPR/-0.94D0, -1.98D0, .1D0,
+ 1 -1.40D0, 0.37D0, .1D0/ ! Re
+ DATA AAPAPI/ 0.3 D0, .267D0,-.01D0,
+ 1 1.66D0, .553D0,-.01D0/ ! Im
+ END
+ SUBROUTINE CKKB ! calculates KK-b scattering amplitude,
+ ! saturated by S*(980) and delta(982) resonances
+ IMPLICIT REAL*8 (A-H,O-Z)
+ COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS,
+ 1 X,Y,Z,T,RP,RPS
+ COMMON/FSI_AAKK/AAKK(9)
+ COMMON/FSI_C/C(10),AM,AMS,DM
+ COMPLEX*16 C
+ S4=AKS+AAKK(1)
+ S=4*S4
+ AKPIPI=DSQRT(S4-AAKK(2))
+ EETA2=(S+AAKK(3)-AAKK(2))**2/4/S
+ AKPIETA=DSQRT(EETA2-AAKK(3))
+ C(1)=AAKK(6)/2/DCMPLX(AAKK(4)-S,
+ ,-AK*AAKK(6)-AKPIPI*AAKK(7))
+ C(1)=C(1)+AAKK(8)/2/DCMPLX(AAKK(5)-S,
+ ,-AK*AAKK(8)-AKPIETA*AAKK(9))
+ RETURN
+ END
+ SUBROUTINE CPAP ! calculates pp-bar scattering amplitude
+ ! accounting for nn-bar->pp-bar channel
+ IMPLICIT REAL*8 (A-H,O-Z)
+ COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
+ COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS,
+ 1 X,Y,Z,T,RP,RPS
+ COMMON/FSI_AAPAP/AAPAPR(3,2),AAPAPI(3,2)
+ COMMON/FSI_C/C(10),AM,AMS,DM
+ COMMON/FSI_2CHA/AK2,AK2S,AAK,HCP2,AMU2_AMU1 ! k* (kappa) for 2-nd channel
+ COMPLEX*16 C
+ DATA AM2/.93956563D0/
+ AMU2_AMU1=AM2/AM ! AM2=2*mu(nn-bar), AM=2*mu(pp-bar)
+ AK2S0=AMU2_AMU1*AKS
+ AK2S =AK2S0-2*AM2*(AM2-AM)
+ IF(AK2S.GE.0.D0)THEN
+ AK2=DCMPLX(DSQRT(AK2S),0.D0) ! k2
+ ELSE
+ AK2=DCMPLX(0.D0,DSQRT(-AK2S)) ! kappa2
+ ENDIF
+ C(10)=C(6+(ISPIN-1)*2)+
+ + DCMPLX(AAPAPR(3,ISPIN)*AKS/2-0.016D0-HCP2,
+ , AAPAPI(3,ISPIN)*AKS/2-AAK) ! (1/f)11
+ C(5)=C(6+(ISPIN-1)*2)+
+ + DCMPLX(AAPAPR(3,ISPIN)*AK2S0/2,
+ , AAPAPI(3,ISPIN)*AK2S0/2)
+ IF(AK2S.GE.0.D0)THEN
+ C(5)=C(5)-DCMPLX(0.D0,AK2)
+ ELSE
+ C(5)=C(5)+DCMPLX(AK2,0.D0) ! (1/f)22
+ ENDIF
+ C(10)=C(10)*C(5)-C(7+(ISPIN-1)*2)*C(7+(ISPIN-1)*2)
+ C(ISPIN)=C(5)/C(10) ! f11
+ C(ISPIN+2)=-C(7+(ISPIN-1)*2)/C(10) ! f12
+ RETURN
+ END
+
+ SUBROUTINE FSIIN(I_ITEST,I_ICH,I_IQS,I_ISI,I_I3C)
+C 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-- total spin S.
+C-- To calculate the CF of two particles (with masses m1, m2 and
+C-- charges C1, C2) the following information is required:
+C-- AM= twice the reduced mass= 2*m1*m2/(m1+m2) in GeV/c^2,
+C-- DM= (m1-m2)/(m1+m2), required if NS=2;
+C-- AC= Bohr radius= 2*137.036*0.1973/(C1*C2*AMH) in fm;
+C-- AC > 1.D9 if C1*C2= 0, AC < 0 if C1*C2 < 0;
+C-- MSPIN= MSPINH(LL)= number of the values of the total pair spin S;
+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-- 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;
+C-- RHO= RHOH(LL,JJ), SF= SFH(LL,JJ), SE= SEH(LL) are spin factors;
+C-- RHO= the probability that the spins j1 and j2 of the two particles
+C-- will combine in a total spin S;
+C-- RHO= (2*S+1)/[(2j1+1)*(2j2+1)] for unpolarized particles;
+C-- RHO= (1-P1*P2)/4 and (3+P1*P2)/4 correspond to S=0 and 1 in the
+C-- case of spin-1/2 particles with polarizations P1 and P2;
+C-----------------------------------------------------------------------
+ IMPLICIT REAL*8 (A-H,O-Z)
+ COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
+ COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS,
+ 1 X,Y,Z,T,RP,RPS
+ COMMON/FSI_SPIN/RHO(10)
+ COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
+ COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
+ 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
+ COMMON/FSI_AA/AA
+ COMMON/FSI_AAPI/AAPI(20,2)/FSI_AAND/AAND(20,4)
+ COMMON/FSI_AAPIN/AAPIN(20,2)
+ COMMON/FSI_AAPAP/AAPAPR(3,2),AAPAPI(3,2)
+ COMMON/FSI_SW/RB(10),EB(10),BK(10),CDK(10),SDK(10),
+ 1 SBKRB(10),SDKK(10)
+ COMMON/FSI_FDH/FDH(30,10),RDH(30,10),EBH(30,10),RBH(30,10)
+ COMMON/FSI_RHOH/RHOH(30,10)
+ COMMON/FSI_AMCH/AM1H(30),AM2H(30),C1H(30),C2H(30),MSPINH(30)
+C============= declarations pour l'appel de READ_FILE()============
+c CHARACTER*10 KEY
+c CHARACTER*8 CH8
+c INTEGER*4 INT4
+c REAL*8 REAL8
+c 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,
+ C 4*.13957D0,4*.493677D0,
+ C 2*1.87561339D0,2*2.80892165D0,2*.497672D0,
+ C 1.87561339D0,3*.93827231D0,.93956563D0,
+ C 1.115684D0,.93827231D0/
+ 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 3*1.115684D0,.93827231D0/
+c--------|---------|---------|---------|---------|---------|---------|----------
+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,2*0.D0, 1.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,3*0.D0,-1.D0/
+C---MSPIN vs (LL)
+ DATA MSPINH/3*2,4*1,2*2,8*1,3,1,2,1,2*1,2*2,1,3*2, 2/
+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,3*.25D0, .25D0,
+ 2 3*.75D0, 4*0.D0, 2*.6667D0, 8*0.D0,
+ 2 .3333D0,.0D0,.75D0,.0D0,2*0.D0,
+ 2 .6667D0,.75D0,0.D0,3*.75D0, .75D0,
+ 3 17*.0D0,.5556D0,3*0.D0, 8*0.D0,1*0.D0,210*0.D0/
+C---Scattering length FD and effective radius RD in fm vs (LL,ISPIN)
+ DATA FDH/
+ 1 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,
+ 1 2.81D0, ! ND potential
+C 1 0.50D0, ! NSC97e potential lam-lam
+ 1 1*0.001D0,
+C c 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.50D0, ! triplet f0 lam-lam=singlet f0 ND
+ ! not contributing in s-wave FSI approx.
+ 2 1*0.001D0,
+ 3 240*0.D0/
+c--------|---------|---------|---------|---------|---------|---------|----------
+ DATA RDH/
+ 1 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,
+ 1 2.95D0, ! ND potential lam-lam
+C 1 10.6D0, ! NSC97e potential lam-lam
+ 1 1*0.D0,
+ 2 3*1.7D0,4*0.D0,2.0D0,2.63D0, 17*0.D0,3.35D0,3.37D0,
+ 2 2.95D0, ! triplet d0 lam-lam=singlet d0 ND
+ ! not contributing in s-wave approx.
+ 2 1*0.D0,
+ 3 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,
+ 1 2.234089D0, ! ND potential lam-lam
+C 1 3.065796D0, ! NSC97e potential lam-lam
+ 1 1*0.001D0,
+ 2 3*2.003144D0,
+ 2 4*0.D0, 2.D0, 4.132163D0, 17*0.D0,
+ 2 2.272703D0, 2.256355D0,
+ 2 2.234089D0, ! triplet potential lam-lam=singlet ND
+ ! not contributing in s-wave FSI approx.
+ 2 1*0.001D0,
+ 3 240*0.D0/
+ DATA EBH/.1149517D0, .1046257D0, .1148757D0, .1186010D0,
+ 1 .7947389D0,2.281208D0,8.7D0,.4D0,.1561219D0,17*0.D0,
+ 1 .1013293D0, .1020966D0,
+ 1 .1080476D0, ! ND potential lam-lam
+C 1 .04115994D0, ! NSC97e potential lam-lam
+ 1 1*0.001D0,
+ 2 3*.1847221D0,
+ 2 4*0.D0, .4D0, .1150687D0, 17*0.D0,
+ 2 .09736083D0, .09708310D0,
+ 2 .1080476D0, ! triplet potential lam-lam= singlet ND
+ ! not contributing in s-wave FSI approx.
+ 2 1*0.001D0,
+ 3 240*0.D0/
+
+C++----- add to be able to call several time-------
+ integer ifirst
+ data ifirst/0/
+ ifirst=ifirst+1
+
+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
+ AC1=1.D10
+ AC2=1.D10
+C=======< condition de calculs >=============
+ NUNIT=11 ! for IBM in Nantes
+C NUNIT=4 ! for SUN in Prague
+C++ CALL readint4(NUNIT,'ITEST ',ITEST)
+C++ CALL readint4(NUNIT,'NS ',NS)
+ ITEST=I_ITEST
+ IF(ITEST.EQ.1)THEN
+C++ CALL readint4(NUNIT,'ICH ',ICH)
+C++ CALL readint4(NUNIT,'IQS ',IQS)
+C++ CALL readint4(NUNIT,'ISI ',ISI)
+C++ CALL readint4(NUNIT,'I3C ',I3C)
+ ICH=I_ICH
+ IQS=I_IQS
+ ISI=I_ISI
+ I3C=I_I3C
+ ENDIF
+C============================================
+ IF (IFIRST.LE.1) THEN
+ DO 3 J1=1,30
+ DO 3 J2=1,10
+ FDH(J1,J2)=FDH(J1,J2)*W
+ RDH(J1,J2)=RDH(J1,J2)*W
+ 3 RBH(J1,J2)=RBH(J1,J2)*W
+C print *,"FD,RD,EB,RB: ",FDH(J1,J2),RDH(J1,J2),EBH(J1,J2),RBH(J1,J2)
+
+ ENDIF
+C===================================
+ RETURN
+ END
+C
+ SUBROUTINE LLINI(lll,I_NS,I_ITEST)
+C===> Initialisation for a given LL value.
+C ===========================================
+ IMPLICIT REAL*8 (A-H,O-Z)
+ COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
+ COMMON/FSI_SPIN/RHO(10)
+ COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
+ COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
+ 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
+ COMMON/FSI_AA/AA
+ COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS,
+ 1 X,Y,Z,T,RP,RPS
+ 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/FSI_FDH/FDH(30,10),RDH(30,10),EBH(30,10),RBH(30,10)
+ COMMON/FSI_RHOH/RHOH(30,10)
+ COMMON/FSI_AMCH/AM1H(30),AM2H(30),C1H(30),C2H(30),MSPINH(30)
+ COMMON/FSI_AAKK/AAKK(9)/FSI_AAPAP/AAPAPR(3,2),AAPAPI(3,2)
+
+C++----- add to be able to call several time-------
+ integer ifirst
+ data ifirst/0/
+ ifirst=ifirst+1
+
+C===> LL - Initialisation ========================================
+C---- Setting particle masses and charges
+ LL=lll
+ NS=I_NS
+ AM1=AM1H(LL)
+ AM2=AM2H(LL)
+ C1=C1H(LL)
+ C2=C2H(LL)
+C print *,"llini: ",AM1,AM2,C1,C2
+C---> Switches:
+C ISI=1(0) the strong interaction between the two particles ON (OFF)
+C IQS=1(0) the quantum statistics ON (OFF);
+C should be OFF for nonidentical particles
+C I3C=1(0) the Coulomb interaction with the nucleus ON (OFF)
+C I3S=1(0) the strong interaction with the nucleus ON (OFF)
+C ICH=1(0) if C1*C2 is different from 0 (is equal to 0)
+C- To switch off the Coulomb force between the two particles
+C put ICH=0 and substitute the strong amplitude parameters by
+C the ones not affected by Coulomb interaction
+C---- --------------------------------------------------------------------
+ IF (I_ITEST.NE.1) 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=1
+ ENDIF
+C---> Calcul. twice the reduced mass (AM), the relative
+C mass difference (DM) and the Bohr radius (AC)
+ AM=2*AM1*AM2/(AM1+AM2)
+ AMS=AM*AM
+ DM=(AM1-AM2)/(AM1+AM2)
+ AC=1.D10
+ C12=C1*C2
+ IF(C12.NE.0.D0)AC=2*137.036D0/(C12*AM)
+C---Setting spin factors
+ MSPIN=MSPINH(LL)
+C print *,"MSPIN: ",MSPIN
+ MSP=MSPIN
+ DO 91 ISPIN=1,10
+ 91 RHO(ISPIN)=RHOH(LL,ISPIN)
+C 91 print *,"RHO: ",ISPIN,RHO(ISPIN)
+C---> Integration limit AA in the spherical wave approximation
+ AA=0.D0
+cc IF(NS.EQ.2.OR.NS.EQ.4)AA=.5D0 !!in 1/GeV --> 0.1 fm
+ IF(NS.EQ.2.OR.NS.EQ.4)AA=6.D0 !!in 1/GeV --> 1.2 fm
+C---> Setting scatt. length (FD), eff. radius (RD) and, if possible,
+C-- also the corresp. square well parameters (EB, RB)
+ DO 55 JJ=1,MSP
+ ISPIN=JJ
+ FD(JJ)=FDH(LL,JJ)
+ RD(JJ)=RDH(LL,JJ)
+ EB(JJ)=EBH(LL,JJ)
+ RB(JJ)=RBH(LL,JJ)
+C print *,"FD,RD,EB,RB: ",FD(JJ),RD(JJ),EB(JJ),RB(JJ)
+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
+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
+C---fm to 1/GeV for pp-bar system
+ IF(LL.EQ.30)THEN
+ IF(IFIRST.LE.1)THEN
+ DO 4 I3=1,3
+ AAPAPR(I3,JJ)=AAPAPR(I3,JJ)*W
+ 4 AAPAPI(I3,JJ)=AAPAPI(I3,JJ)*W
+C 4 print *,"AAPAPR,AAPAPI: ",AAPAPR(I3,JJ),AAPAPI(I3,JJ)
+C--- Calculates complex elements M11=M22=C(6), M12=M21=C(7) for I=0
+C--- at k*=0 M11=M22=C(8), M12=M21=C(9) for I=1
+ C(7+(JJ-1)*2)=2*DCMPLX(AAPAPR(1,JJ),AAPAPI(1,JJ))*
+ * DCMPLX(AAPAPR(2,JJ),AAPAPI(2,JJ)) ! 2a_0Sa_1S
+ C(6+(JJ-1)*2)=DCMPLX(AAPAPR(1,JJ)+AAPAPR(2,JJ),
+ , AAPAPI(1,JJ)+AAPAPI(2,JJ))/
+ / C(7+(JJ-1)*2) ! M11=M22
+ C(7+(JJ-1)*2)=-DCMPLX(AAPAPR(1,JJ)-AAPAPR(2,JJ),
+ , AAPAPI(1,JJ)-AAPAPI(2,JJ))/
+ / C(7+(JJ-1)*2) ! M12=M21
+ ENDIF
+ ENDIF
+C---Calculation continues for any system (any LL)
+ 55 CONTINUE
+ RETURN
+ END
+C=======================================================
+C
+
+c++ This routine is used to init mass and charge of the nucleus.
+
+ SUBROUTINE FSINUCL(R_AMN,R_CN)
+
+ IMPLICIT REAL*8 (A-H,O-Z)
+ COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
+
+ AMN=R_AMN
+ CN=R_CN
+
+ RETURN
+ END
+
+C======================================================
+C
+
+ SUBROUTINE FSIMOMENTUM(PP1,PP2)
+
+ IMPLICIT REAL*8 (A-H,O-Z)
+ COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1, ! particle momenta in NRF
+ 1 P2X,P2Y,P2Z,E2,P2
+
+
+ REAL*8 PP1(3),PP2(3)
+c Print *,"momentum",pp1,pp2
+ P1X=PP1(1)
+ P1Y=PP1(2)
+ P1Z=PP1(3)
+ P2X=PP2(1)
+ P2Y=PP2(2)
+ P2Z=PP2(3)
+ RETURN
+ END
+
+
+
+C======================================================
+C
+
+ SUBROUTINE FSIPOSITION(XT1,XT2)
+
+ IMPLICIT REAL*8 (A-H,O-Z)
+ COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, !4-coord. of emis. points in NRF
+ 1 X2,Y2,Z2,T2,R2
+
+ REAL*8 XT1(4),XT2(4)
+clc print *,'fsi',xt1,xt2
+ X1=XT1(1)
+ Y1=XT1(2)
+ Z1=XT1(3)
+ T1=XT1(4)
+ X2=XT2(1)
+ Y2=XT2(2)
+ Z2=XT2(3)
+ T2=XT2(4)
+ RETURN
+ END
+
+
+C======================================================
+C======================================================
+C
+ subroutine BoostToPrf()
+ 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_P12/P12X,P12Y,P12Z,E12,P12,AM12,EPM
+
+ 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)
+CW WRITE(6,38)'RP ',RP,'X ',X,Y,Z,T
+38 FORMAT(A7,E11.4,A7,4E11.4)
+
+ CVK=(P12X*PPX+P12Y*PPY+P12Z*PPZ)/(P12*AK)
+ V=P12/E12
+ return
+ end
+
+ 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/FSI_2CHA/AK2,AK2S,AAK,HCP2,AMU2_AMU1 ! k* (kappa) for 2-nd channel
+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)
+CW WRITE(6,38)'RP ',RP,'X ',X,Y,Z,T
+38 FORMAT(A7,E11.4,A7,4E11.4)
+
+ CVK=(P12X*PPX+P12Y*PPY+P12Z*PPZ)/(P12*AK)
+ V=P12/E12
+
+C ACH=1.D0
+ 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
+C AAK=ACH*AK !
+C HCP2=2*HCP/AC ! needed to calculate C(JJ) for charged particles
+ 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=6-26 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
+ IF(LL.EQ.30)GOTO 28 ! pp-bar
+ 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
+c IF(LL.EQ.17)C(JJ)=DCMPLX(-3.29D0,3.55D0) ! Martin'76 (K-p)
+c IF(LL.EQ.17)C(JJ)=DCMPLX(3.29D0,3.55D0) ! Martin'76 (K-p) WRONG SIGN!!!
+ IF(LL.EQ.17)C(JJ)=DCMPLX(-2.585D0,4.156D0) ! Borasoy'04 (K-p)
+c IF(LL.EQ.17)C(JJ)=DCMPLX(-3.371D0,3.244D0) ! Martin'81 (K-p)
+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 ! needed to calculate C(JJ) for charged particles
+ C(JJ)=1/(1/C(JJ)-HCP2+DCMPLX(0.D0,AK-AAK))
+ GOTO 30
+ 28 CONTINUE
+C---Calc. pp-bar s-wave scatt. ampl.
+ CALL CPAP
+ 30 CONTINUE
+C***********************************************************************
+ CALL VZ(WEI)
+ RETURN
+ END
+ SUBROUTINE VZ(WEI)
+C==> Calculates the weight WEI due to FSI
+ IMPLICIT REAL*8 (A-H,O-Z)
+ COMMON/FSI_JR/JRAT
+ COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS,
+ 1 X,Y,Z,T,RP,RPS
+ COMMON/FSI_SPIN/RHO(10)
+ COMMON/FSI_ETA/ETA
+ COMMON/FSI_AA/AA
+ COMMON/FSI_FFF/F12,F21
+ COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
+ COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
+ COMMON/FSI_FD/FD(10),RD(10)
+ COMMON/FSI_RR/F(10)
+ COMMON/FSI_C/C(10),AM,AMS,DM
+ COMMON/FSI_COULPH/EIDC
+ 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
+ COMMON/FSI_2CHA/AK2,AK2S,AAK,HCP2,AMU2_AMU1 ! k* (kappa) for 2-nd channel
+ WEI=0.D0
+ IF(JRAT.EQ.1)GOTO 11
+ RHOS=AK*RP
+ HS=X*PPX+Y*PPY+Z*PPZ
+ 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
+ Z8=CMPLX(1.,SNGL(ETA))
+ Z8=CGAMMA(Z8)
+ EIDC=Z8/CABS(Z8)
+ 2 CALL FF(RHOS,HS)
+ 11 MSP=MSPIN
+ 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
+ IF(IQS.EQ.0)GOTO 5 ! the quantum statistics ON (OFF) if IQS=1(0)
+ JSIGN=-1
+ DO 1 JJ=1,MSP
+ JSIGN=-JSIGN
+ G=F(JJ)*C(JJ)
+ IF(ICH.EQ.1)G=G*ACHR
+ PSI12=FF12*(F12+G)
+ PSI21=FF21*(F21+G)
+ G=PSI12+JSIGN*PSI21
+ 1 WEI=WEI+RHO(JJ)*(DREAL(G)**2+DIMAG(G)**2)
+ GOTO 8
+ 5 DO 6 JJ=1,MSP
+ G=F(JJ)*C(JJ)
+ IF(ICH.EQ.1)G=G*ACHR
+CW WRITE(6,38)'JJ ',JJ,'F ',DREAL(F(JJ)),DIMAG(F(JJ))
+CW WRITE(6,38)'JJ ',JJ,'C ',DREAL(C(JJ)),DIMAG(C(JJ))
+CW WRITE(6,38)'JJ ',JJ,'G ',DREAL(G),DIMAG(G)
+CW WRITE(6,38)'JJ ',JJ,'F12+G ',DREAL(F12+G),DIMAG(F12+G)
+CW WRITE(6,38)'JJ ',JJ,'F21+G ',DREAL(F21+G),DIMAG(F21+G)
+38 FORMAT(A7,I3,A7,2E11.4)
+ PSI12=FF12*(F12+G)
+ 6 WEI=WEI+RHO(JJ)*(DREAL(PSI12)**2+DIMAG(PSI12)**2)
+c--- Account for nn-bar->pp-bar channel ---------------------------
+ IF(LL.EQ.30)THEN
+ DO 61 JJ=1,MSP
+ HH=RHO(JJ)*(DREAL(C(JJ+2))**2+DIMAG(C(JJ+2))**2)*
+ * AMU2_AMU1*ACH/RPS
+ IF(AK2S.LT.0)HH=HH*DEXP(-2*RP*AK2)
+ 61 WEI=WEI+HH
+ ENDIF
+c------------------------------------------------------------------
+ RETURN
+ 4 PSI12=FF12*F12
+ IF(IQS.EQ.0)GOTO 50 ! the quantum statistics ON (OFF) if IQS=1(0)
+ PSI21=FF21*F21
+ JSIGN=-1
+ DO 3 JJ=1,MSP
+ JSIGN=-JSIGN
+ G=PSI12+JSIGN*PSI21
+ 3 WEI=WEI+RHO(JJ)*(DREAL(G)**2+DIMAG(G)**2)
+ GOTO 8
+ 50 WEI=DREAL(PSI12)**2+DIMAG(PSI12)**2
+ RETURN
+ 8 WEI=WEI/2
+ RETURN
+ END
+ SUBROUTINE FIRT
+C---CALC. THE F(JJ)
+C-- F(JJ)*C(JJ)= DEVIATION OF THE BETHE-SALPETER AMPL. FROM PLANE WAVE
+ IMPLICIT REAL*8 (A-H,O-Z)
+ COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS,
+ 1 X,Y,Z,T,RP,RPS
+ COMMON/FSI_SHH/SH,CHH
+ COMMON/FSI_BP/B,P
+ COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
+ COMMON/FSI_C/C(10),AM,AMS,DM
+ 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)
+ EQUIVALENCE(RSS,RP),(TSS,T)
+ COMPLEX*16 F,C,CH1
+ MSP=MSPIN
+ DO 10 JJ=1,MSP
+ IF(JJ.GT.1)GOTO 3
+ XRA=2*RSS/AC
+ IF(AK.NE.0.D0)GOTO2
+ SHK=1.D0
+ SH=.0D0
+ SHH=SH
+ CHH=1/RSS
+ GOTO3
+ 2 H=AK*RSS
+ CALL GST(XRA,H)
+ SH=SH/RSS
+ CHH=CHH/RSS
+ SHH=SH
+ IF(ICH.EQ.1) SHH=ACH*SH
+ 3 IF(NS.EQ.2)GOTO1
+C---F= ASYMPTOTIC FORMULA (T= 0 APPROX.); NS= 4
+ 6 F(JJ)=DCMPLX(CHH,SHH)
+ IF(NS.NE.1)GOTO 10
+C---F INSIDE THE SQUARE-WELL (T= 0 APPROX.); NS= 1
+ IF(RSS.GE.RB(JJ)) GOTO 10
+ IF(AK.NE.0.D0.AND.JJ.EQ.1)SHK=SH/AK
+ H=BK(JJ)*RSS
+ CALL GST(XRA,H)
+ SKR=B*BK(JJ)
+ F(JJ)=DCMPLX(CDK(JJ),SDK(JJ))*SKR
+ CH1=(SDKK(JJ)*SKR-SHK*SBKRB(JJ))/C(JJ)
+ F(JJ)=(F(JJ)+CH1)/SBKRB(JJ)
+ GOTO 10
+ 1 CONTINUE
+C---F= ASYMPTOTIC FORMULA (T= 0 NOT REQUIRED); NS= 2
+ IF(JJ.GT.1)GOTO 8
+ IF(TSS.EQ.0.D0)GOTO6
+ TSSA=DABS(TSS)
+ IF(DM.NE.0.D0)GOTO 11
+ H=AM*.5D0/TSSA
+ 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
+ 4 CONTINUE
+ H1=AK*TSSA/AM
+ 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
+ FS2=0.D0
+ FC1=0.D0
+ FC2=0.D0
+ DO 13 I=1,2
+ IF(I.EQ.1)TSSH=TSSA*(1+DM)
+ IF(I.EQ.2)TSSH=TSSA*(1-DM)
+ H=AM*.5D0/TSSH
+ 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
+ GOTO 13
+ 12 CONTINUE
+ H1=AK*TSSH/AM
+ 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
+ A12=FS1-FC2
+ A21=FS2-FC1
+ A2=.5D0*(CHH*(A12+A21)+SH*(A12-A21))+SHH
+ A1=.5D0*(CHH*(C12+S12)+SH*(C12-S12))
+ F(JJ)=.3989422D0*DCMPLX(A1,A2)
+ GOTO 10
+ 8 F(JJ)=F(1)
+ 10 CONTINUE
+ RETURN
+ END
+ FUNCTION EXF(X)
+ IMPLICIT REAL*8 (A-H,O-Z)
+ IF(X.LT.-15.D0) GO TO 1
+ EXF=DEXP(X)
+ RETURN
+ 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.
+ IMPLICIT REAL*8 (A-H,O-Z)
+ COMMON/FSI_BP/B,P
+ DIMENSION BH(3),PH(3)
+ DATA ERR/1.D-7/
+ BH(1)=1.D0
+ PH(1)=1.D0
+ PH(2)=.0D0
+ BH(2)=.5D0*X
+ B=1+BH(2)
+ P=1.D0
+ HS=H*H
+ J=0
+ 2 J=J+1
+ BH(3)=(X*BH(2)-HS*BH(1))/((J+1)*(J+2))
+ PH(3)=(X*PH(2)-HS*PH(1)-(2*J+1)*X*BH(2))/(J*(J+1))
+ B=B+BH(3)
+ P=P+PH(3)
+ Z=DABS(BH(2))+DABS(BH(3))+DABS(PH(2))+DABS(PH(3))
+ IF(Z.LT.ERR)RETURN
+ BH(1)=BH(2)
+ BH(2)=BH(3)
+ PH(1)=PH(2)
+ 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.
+ IMPLICIT REAL*8 (A-H,O-Z)
+ COMMON/FSI_BP/B,P
+ COMMON/FSI_SHH/SH,CHH
+ COMMON/FSI_ETA/ETA
+ COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
+ COMMON/FSI_COULPH/EIDC
+ COMPLEX*16 EIDC,GST
+ ARG=H-ETA*DLOG(2*H)
+ GST=DCMPLX(DCOS(ARG),DSIN(ARG))
+ GST=ACHR*EIDC*GST
+ CHH=DREAL(GST)
+ SH=DIMAG(GST)/ACH
+ B=SH/H
+ RETURN
+ END
+ SUBROUTINE FF(RHO,H)
+C---Calc. F12, F21;
+C-- F12= FF0* plane wave, FF0=F*ACHR,
+C---F is the confluent hypergeometric function,
+C-- ACHR=sqrt(ACH), where ACH is the Coulomb factor
+ IMPLICIT REAL*8 (A-H,O-Z)
+ COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
+ COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
+ COMMON/FSI_ETA/ETA
+ COMMON/FSI_FFF/F12,F21
+ COMPLEX*16 FF0,F12,F21
+ C=DCOS(H)
+ S=DSIN(H)
+ F12=DCMPLX(C,-S)
+ F21=DCMPLX(C,S)
+ IF(ICH.EQ.0)RETURN
+ RHOP=RHO+H
+ RHOM=RHO-H
+ F12=FF0(RHO,H)*F12
+ 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
+C-- ACHR=sqrt(ACH), where ACH is the Coulomb factor
+ IMPLICIT REAL*8 (A-H,O-Z)
+ COMPLEX*16 FAS,EIDC,ZZ1
+ COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
+ COMMON/FSI_ETA/ETA
+ COMMON/FSI_COULPH/EIDC
+ D1=DLOG(RKS)*ETA
+ D2=ETA*ETA/RKS
+ ZZ1=DCMPLX(DCOS(D1),DSIN(D1))/EIDC
+ FAS=DCMPLX(1.D0,-D2)*ZZ1
+ 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
+C-- (Eq. (15) of G-K-L-L), F= 1 at r* << AC
+C-- ACHR=sqrt(ACH), where ACH is the Coulomb factor
+ IMPLICIT REAL*8 (A-H,O-Z)
+ COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
+ COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
+ COMMON/FSI_ETA/ETA
+ COMPLEX*16 ALF,ALF1,Z,S,A,FF0,FAS
+ DATA ERR/1.D-5/
+ S=DCMPLX(1.D0,0.D0)
+ FF0=S
+ RHOP=RHO+H
+CC GOTO 5 ! rejects the approx. calcul. of hyperg. f-ion F
+ IF(RHOP.LT.20.D0)GOTO5
+ FF0=FAS(RHOP) ! approx. calc.
+ RETURN
+ 5 ALF=DCMPLX(.0D0,-ETA)
+ ALF1=ALF-1
+ Z=DCMPLX(.0D0,RHOP)
+ J=0
+ 3 J=J+1
+ A=(ALF1+J)/(J*J)
+ S=S*A*Z
+ FF0=FF0+S
+ ZR=DABS(DREAL(S))
+ ZI=DABS(DIMAG(S))
+ IF((ZR+ZI).GT.ERR)GOTO3
+ 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
+C-- the gamma function)
+ IMPLICIT REAL*8 (A-H,O-Z)
+ DIMENSION BN(15)
+ DATA BN/.8333333333D-1,.8333333333D-2,.396825396825D-2,
+ 1 .4166666667D-2,.7575757576D-2,.2109279609D-1,
+ 2 .8333333333D-1,.4432598039D0 ,.305395433D1,
+ 3 .2645621212D2, .2814601449D3, .3607510546D4,
+ 4 .5482758333D5, .9749368235D6, .200526958D8/
+ X=DABS(XA)
+ IF(X.LT..33D0) GOTO 1
+CC IF(X.GE.3.5D0) GO TO 2
+ S=.0D0
+ N=0
+ 3 N=N+1
+ DS=1.D0/N/((N*X)**2+1)
+ S=S+DS
+ IF(DS.GT.0.1D-12) GOTO 3
+C---Provides 7 digit accuracy
+ HC=S-.5772156649D0+DLOG(X)
+ RETURN
+CC 2 HC=1.2D0/X**2+DLOG(X)-.5772156649 D0
+CC RETURN
+ 1 X2=X*X
+ XP=X2
+ HC=0.D0
+ IMA=9
+ IF(X.LT.0.1D0)IMA=3
+ DO 4 I=1,IMA
+ HC=HC+XP*BN(I)
+ 4 XP=X2*XP
+ RETURN
+ END
+ FUNCTION ACP(X)
+C--- ACP = COULOMB PENETRATION FACTOR
+ IMPLICIT REAL*8 (A-H,O-Z)
+ IF(X.LT.0.05D0.AND.X.GE.0.D0) GO TO 1
+ Y=6.2831853D0/X
+ ACP=Y/(EXF(Y)-1)
+ RETURN
+ 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).
+ IMPLICIT REAL*8 (A-H,O-Z)
+ COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
+ COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
+ COMMON/FSI_SHH/SH,CHH
+ COMMON/FSI_BP/B,P
+ 1 IF(ICH.EQ.1)GOTO 2
+ 3 SH=DSIN(H)
+ CHH=DCOS(H)
+ B=1.D0
+ IF(H.NE.0.D0)B=SH/H
+ P=CHH
+ RETURN
+ 2 CONTINUE
+ IF(H.GT.15.D0)GOTO4 ! comment out if you want to reject
+ ! the approximate calculation of hyperg. f-ion G
+ CALL SEQ(X,H) ! exact calculation
+ SH=H*B
+ CHH=P+B*X*(DLOG(DABS(X))+HPR)
+ RETURN
+ 4 CALL SEQA(X,H)
+ RETURN
+ END
+ FUNCTION FF1(RHO,H)
+C---FF1=FF0; used for particle-nucleus system
+C-- FF0=F12*ACHR
+C-- F12 is the confluent hypergeometric function
+C-- (Eq. (15) of G-K-L-L), F12= 1 at r* << AC
+C-- ACHR=sqrt(ACH), where ACH is the Coulomb factor
+ IMPLICIT REAL*8 (A-H,O-Z)
+ COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
+ COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
+ COMMON/FSI_ETA/ETA
+ COMMON/FSI_COULPH/EIDC
+ COMMON/FSI_ICH1/ICH1
+ COMPLEX*16 FF0,FF1
+ COMPLEX*16 EIDC
+ COMPLEX*8 Z8,CGAMMA
+ FF1=DCMPLX(1.D0,0.D0)
+ IF(ICH1.EQ.0)GOTO 2
+ 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),
+ 1 SBKRB(10),SDKK(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
+ ACH=1.D0
+ IF(ICH.EQ.0)GOTO 1
+ XH=AC*AK
+ HCP=HC(XH)
+ HPR=HCP+.1544313298D0
+ ACH=ACP(XH)
+ HAC=2*HCP/AC
+ 1 AKS=AK**2
+ BK(JJ)=DSQRT(AKS+EB(JJ)**2) ! kappa=kp
+ X=2*RB(JJ)/AC
+ H=BK(JJ)*RB(JJ) ! kp*d
+ CALL GST(X,H)
+ BRHO=B ! B(kp,d)
+ SBKRB(JJ)=SH ! kp*d*B(kp,d)
+ CALL DERIW(X,H)
+ BRHOP=BPR ! B'(kp,d)= dB(kp,r)/dln(r) at r=d
+ H=AK*RB(JJ)
+ CALL GST(X,H)
+ CDK(JJ)=CHH ! ReG(k,d)
+ BRHOS=B ! B(k,d)
+ PRHOS=P ! P(k,d)
+ SDK(JJ)=SH
+ IF(ICH.EQ.0)GOTO 2
+ SDK(JJ)=ACH*SH ! ImG(k,d)
+ IF(AK.EQ.0.D0.AND.AC.LT.0.D0)SDK(JJ)=3.14159*X*B
+ 2 SDKK(JJ)=RB(JJ)
+ IF(AK.NE.0.D0)SDKK(JJ)=SH/AK ! d*B(k,d)
+ CALL DERIW(X,H) ! PPR=P'(k,d)= dP(k,r)/dln(r) at r=d
+ ZZ=PPR-PRHOS
+ IF(ICH.EQ.1)ZZ=ZZ+X*(BRHOS+BPR*(DLOG(DABS(X))+HPR))
+C ZZ= P'(k,d)-P(k,d)+x*{B(k,d)+B'(k,d)*[ln!x!+2*C-1+h(k*ac)]}
+ GG=(BRHOP*CDK(JJ)-BRHO*ZZ)/RB(JJ)
+C GG= [B'(kp,d)*ReG(k,d)-B(kp,d)*ZZ]/d
+ G=GG/(BRHO*BPR-BRHOP*BRHOS)
+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)
+ COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
+ COMMON/FSI_BP/B,P/FSI_DERIV/BPR,PPR
+ HH=.1D-3
+ CALL GST(X,H-HH)
+ Q1=P
+ B1=B
+ CALL GST(X,H+HH)
+ HHH=HH+HH
+ BPR=H*(B-B1)/HHH
+ PPR=H*(P-Q1)/HHH
+ IF(ICH.EQ.0)RETURN
+ CALL GST(X-HH,H)
+ Q1=P
+ B1=B
+ CALL GST(X+HH,H)
+ BPR=BPR+X*(B-B1)/HHH
+ PPR=PPR+X*(P-Q1)/HHH
+ RETURN
+ END
+C================================================================
+
+ SUBROUTINE READ_FILE(KEY,CH8,INT4,REAL8,IERR,NUNIT)
+C ==========
+C
+C Routine to read one parameter of the program in the file
+C DATA NUNIT defined in FSI3B EXEC
+C NUNIT=11 for IBM in Nantes, 4 for SUN in Prague
+C
+C INPUT : KEY (CHARACTER*10) :
+C OUTPUT : case of KEY : CH8 : (CHARACTER*8)
+C INT4 : (INTEGER*4)
+C REAL8 : (REAL*8)
+C (only one of them)
+C IERR (INTEGER) : 0 : no error
+C 1 : key not found
+
+ CHARACTER*10 KEY,TEST
+ CHARACTER*4 TYPE
+ CHARACTER*8 CH8
+ INTEGER*4 INT4
+ REAL*8 REAL8
+ INTEGER*4 IERR
+
+ IERR=0
+ REWIND(NUNIT)
+1 READ(NUNIT,FMT='(A10,2X,A4)')TEST,TYPE
+ IF (TEST.EQ.KEY) THEN
+ BACKSPACE(NUNIT)
+ IF (TYPE.EQ.'CHAR') READ(NUNIT,FMT='(18X,A8,54X)')CH8
+ IF (TYPE.EQ.'INT4') READ(NUNIT,FMT='(18X,I8,54X)')INT4
+ IF (TYPE.EQ.'REA8') READ(NUNIT,FMT='(18X,F10.5,52X)')REAL8
+ ELSE
+ IF (TEST.NE.'* E.O.F. *') THEN
+ GOTO 1
+ ELSE
+ IERR=1
+ ENDIF
+ ENDIF
+c IF(IERR.EQ.1)STOP
+ RETURN
+ END
+
--- /dev/null
+///////////////////////////////////////////////////////////////////////////
+// //
+// AliFemtoModelWeightGeneratorLednicky : the most advanced weight //
+// generator available. Supports a large number of different pair types //
+// and interaction types. Can calculate pair weights coming from //
+// quantum statistics, coulomb interation and strong interaction ot any //
+// combination of the three, as applicable. //
+// This class is a wrapper for the fortran code provided by Richard //
+// Lednicky. //
+// //
+///////////////////////////////////////////////////////////////////////////
+
+//#include "StHbtMaker/ThCorrFctn/AliFemtoModelWeightGeneratorLednicky.h"
+#include "AliFemtoModelWeightGeneratorLednicky.h"
+#include "AliFemtoModelHiddenInfo.h"
+//#include "StarCallf77.h"
+//#include <strstream.h>
+//#include <iomanip.h>
+//#include <stream>
+//#include <iomanip>
+#include <sstream>
+
+#ifdef SOLARIS
+# ifndef false
+typedef int bool;
+#define false 0
+#define true 1
+# endif
+#endif
+
+#ifdef WIN32
+# ifdef CERNLIB_MSSTDCALL
+# define F77_UCASE
+# define type_of_call _stdcall
+# ifndef CERNLIB_QXCAPT
+# define CERNLIB_QXCAPT
+# endif
+# else
+# define F77_LCASE
+# ifndef CERNLIB_QXNO_SC
+# define CERNLIB_QXNO_SC
+# endif
+# endif
+# define type_of_call _stdcall
+# define DEFCHARD const char* , const int
+# define DEFCHARL
+# define PASSCHARD(string) string, strlen(string)
+# define PASSCHARL(string)
+#else
+# define DEFCHARD const char*
+# define DEFCHARL , const int
+# define PASSCHARD(string) string
+# define PASSCHARL(string) , strlen(string)
+#endif
+#ifdef CERNLIB_QXCAPT
+# define F77_NAME(name,NAME) NAME
+#else
+# if defined(CERNLIB_QXNO_SC)
+# define F77_NAME(name,NAME) name
+# else
+# define F77_NAME(name,NAME) name##_
+# endif
+#endif
+#ifndef type_of_call
+# define type_of_call
+#endif
+
+// --- Prototype of the function used in the weight calculator
+// (in FsiWeightLedinicky.F)
+#define fsiin F77_NAME(fsiin,FSIIN)
+extern "C" {void type_of_call F77_NAME(fsiin,FSIIN)(const int &itest,const int &ich, const int &iqs, const int &isi,const int &i3c);}
+#define llini F77_NAME(llini,LLINI)
+extern "C" {void type_of_call F77_NAME(llini,LLINI)(const int &lll,const int &ns, const int &itest);}
+
+#define fsinucl F77_NAME(fsinucl,FSINUCL)
+extern "C" {void type_of_call F77_NAME(fsinucl,FSINUCL)(const double &mn,const double &cn);}
+#define fsimomentum F77_NAME(fsimomentum,FSIMOMENTUM)
+extern "C" {void type_of_call F77_NAME(fsimomentum,FSIMOMENTUM)(double &p1,double &p2);}
+#define fsiposition F77_NAME(fsiposition,FSIPOSITION)
+extern "C" {void type_of_call F77_NAME(fsiposition,FSIPOSITION)(double &x1,double &x2);}
+#define fsiw F77_NAME(fsiw,FSIW)
+extern "C" {void type_of_call F77_NAME(fsiw,FSIW)(const int &i,double &weif,
+ double &wei,double &wein);}
+#define ltran12 F77_NAME(ltran12,LTRAN12)
+extern "C" {void type_of_call ltran12_();}
+
+// Test function for Lambda potential
+//#define printlam F77_NAME(printlam,PRINTLAM)
+//extern "C" {void type_of_call printlam_();}
+//there is not PRINTLAM in *.F file
+
+// --- Additional prototyping of some CERN functions (in FsiTool.F)
+typedef float REAL;
+typedef struct { REAL re; REAL im; } COMPLEX;
+#define cgamma F77_NAME(cgamma,CGAMMA)
+extern "C" {COMPLEX type_of_call cgamma_(COMPLEX*);}
+
+#ifdef __ROOT__
+ClassImp(AliFemtoModelWeightGeneratorLednicky)
+#endif
+
+AliFemtoModelWeightGeneratorLednicky::AliFemtoModelWeightGeneratorLednicky() :
+ AliFemtoModelWeightGenerator(),
+ fWei(0), fWein(0), fWeif(0), fWeightDen(0),
+ fItest(0),fIch(1),fIqs(1),fIsi(1),fI3c(0),
+ fNuclMass(1.),fNuclCharge(0.),
+ fSphereApp(false),fT0App(false) ,
+ fLL(0), fNuclChargeSign(1), fSwap(0), fLLMax(30), fLLName(0),
+ fNumProcessPair(0), fNumbNonId(0)
+{
+ // default constructor
+ fLLName=new char*[fLLMax+1];
+ fNumProcessPair=new int[fLLMax+1];
+ int i;
+ for (i=1;i<=fLLMax;i++) {fLLName[i]=new char[40];fNumProcessPair[i]=0;}
+ strcpy( fLLName[1],"neutron neutron");
+ strcpy( fLLName[2],"proton proton");
+ strcpy( fLLName[3],"neutron proton");
+ strcpy( fLLName[4],"alpha alpha");
+ strcpy( fLLName[5],"pi+ pi-");
+ strcpy( fLLName[6],"pi0 pi0");
+ strcpy( fLLName[7],"pi+ pi+");
+ strcpy( fLLName[8],"neutron deuteron");
+ strcpy( fLLName[9],"proton deuteron");
+ strcpy( fLLName[10],"pi+ K-");
+ strcpy( fLLName[11],"pi+ K+");
+ strcpy( fLLName[12],"pi+ proton");
+ strcpy( fLLName[13],"pi- proton");
+ strcpy( fLLName[14],"K+ K-");
+ strcpy( fLLName[15],"K+ K+");
+ strcpy( fLLName[16],"K+ proton");
+ strcpy( fLLName[17],"K- proton");
+ strcpy( fLLName[18],"deuteron deuteron");
+ strcpy( fLLName[19],"deuton alpha");
+ strcpy( fLLName[20],"triton triton");
+ strcpy( fLLName[21],"triton alpha");
+ strcpy( fLLName[22],"K0 K0");
+ strcpy( fLLName[23],"K0 K0b");
+ strcpy( fLLName[24],"deuteron triton");
+ strcpy( fLLName[25],"proton triton");
+ strcpy( fLLName[26],"proton alpha");
+ strcpy( fLLName[27],"proton lambda");
+ strcpy( fLLName[28],"neutron lambda");
+ strcpy( fLLName[29],"Lambda lambda");// gael 21May02
+ strcpy( fLLName[30],"Proton Anti-proton");// gael 21May02
+ FsiInit();
+ FsiNucl();
+};
+//______________________
+AliFemtoModelWeightGeneratorLednicky::AliFemtoModelWeightGeneratorLednicky(const AliFemtoModelWeightGeneratorLednicky &aWeight):
+ AliFemtoModelWeightGenerator(),
+ fWei(0), fWein(0), fWeif(0), fWeightDen(0),
+ fItest(0),fIch(1),fIqs(1),fIsi(1),fI3c(0),
+ fNuclMass(1.),fNuclCharge(0.),
+ fSphereApp(false),fT0App(false) ,
+ fLL(0), fNuclChargeSign(1), fSwap(0), fLLMax(30), fLLName(0),
+ fNumProcessPair(0), fNumbNonId(0)
+{
+ // copy constructor
+ fWei = aWeight.fWei;
+ fWein = aWeight. fWein;
+ fWeif = aWeight. fWeif;
+ fWeightDen = aWeight.fWeightDen;
+
+ fItest = aWeight.fItest;
+ fIch = aWeight.fIch;
+ fIqs = aWeight.fIqs;
+ fIsi = aWeight.fIsi;
+ fI3c = aWeight.fI3c;
+ fNuclMass = aWeight.fNuclMass;
+ fNuclCharge = aWeight.fNuclCharge;
+ fSphereApp = aWeight.fSphereApp;
+ fT0App = aWeight.fT0App;
+ fLL = aWeight.fLL;
+ fNuclChargeSign = aWeight.fNuclChargeSign;
+ fSwap = aWeight.fSwap;
+ fLLName = aWeight.fLLName;
+ fNumProcessPair = aWeight.fNumProcessPair;
+ fNumbNonId = aWeight.fNumbNonId;
+ fLLName=new char*[fLLMax+1];
+ fNumProcessPair=new int[fLLMax+1];
+ int i;
+ for (i=1;i<=fLLMax;i++) {fLLName[i]=new char[40];fNumProcessPair[i]=0;}
+ strcpy( fLLName[1],"neutron neutron");
+ strcpy( fLLName[2],"proton proton");
+ strcpy( fLLName[3],"neutron proton");
+ strcpy( fLLName[4],"alpha alpha");
+ strcpy( fLLName[5],"pi+ pi-");
+ strcpy( fLLName[6],"pi0 pi0");
+ strcpy( fLLName[7],"pi+ pi+");
+ strcpy( fLLName[8],"neutron deuteron");
+ strcpy( fLLName[9],"proton deuteron");
+ strcpy( fLLName[10],"pi+ K-");
+ strcpy( fLLName[11],"pi+ K+");
+ strcpy( fLLName[12],"pi+ proton");
+ strcpy( fLLName[13],"pi- proton");
+ strcpy( fLLName[14],"K+ K-");
+ strcpy( fLLName[15],"K+ K+");
+ strcpy( fLLName[16],"K+ proton");
+ strcpy( fLLName[17],"K- proton");
+ strcpy( fLLName[18],"deuteron deuteron");
+ strcpy( fLLName[19],"deuton alpha");
+ strcpy( fLLName[20],"triton triton");
+ strcpy( fLLName[21],"triton alpha");
+ strcpy( fLLName[22],"K0 K0");
+ strcpy( fLLName[23],"K0 K0b");
+ strcpy( fLLName[24],"deuteron triton");
+ strcpy( fLLName[25],"proton triton");
+ strcpy( fLLName[26],"proton alpha");
+ strcpy( fLLName[27],"proton lambda");
+ strcpy( fLLName[28],"neutron lambda");
+ strcpy( fLLName[29],"Lambda lambda");// gael 21May02
+ strcpy( fLLName[30],"Proton Anti-proton");// gael 21May02
+ FsiInit();
+ FsiNucl();
+}
+
+AliFemtoModelWeightGeneratorLednicky& AliFemtoModelWeightGeneratorLednicky::operator=(const AliFemtoModelWeightGeneratorLednicky& aWeight)
+{
+ // assignment operator
+ if (this == &aWeight)
+ return *this;
+
+ fWei = aWeight.fWei;
+ fWein = aWeight. fWein;
+ fWeif = aWeight. fWeif;
+ fWeightDen = aWeight.fWeightDen;
+
+ fItest = aWeight.fItest;
+ fIch = aWeight.fIch;
+ fIqs = aWeight.fIqs;
+ fIsi = aWeight.fIsi;
+ fI3c = aWeight.fI3c;
+ fNuclMass = aWeight.fNuclMass;
+ fNuclCharge = aWeight.fNuclCharge;
+ fSphereApp = aWeight.fSphereApp;
+ fT0App = aWeight.fT0App;
+ fLL = aWeight.fLL;
+ fNuclChargeSign = aWeight.fNuclChargeSign;
+ fSwap = aWeight.fSwap;
+ fLLName = aWeight.fLLName;
+ fNumProcessPair = aWeight.fNumProcessPair;
+ fNumbNonId = aWeight.fNumbNonId;
+ fLLName=new char*[fLLMax+1];
+ fNumProcessPair=new int[fLLMax+1];
+ int i;
+ for (i=1;i<=fLLMax;i++) {fLLName[i]=new char[40];fNumProcessPair[i]=0;}
+ strcpy( fLLName[1],"neutron neutron");
+ strcpy( fLLName[2],"proton proton");
+ strcpy( fLLName[3],"neutron proton");
+ strcpy( fLLName[4],"alpha alpha");
+ strcpy( fLLName[5],"pi+ pi-");
+ strcpy( fLLName[6],"pi0 pi0");
+ strcpy( fLLName[7],"pi+ pi+");
+ strcpy( fLLName[8],"neutron deuteron");
+ strcpy( fLLName[9],"proton deuteron");
+ strcpy( fLLName[10],"pi+ K-");
+ strcpy( fLLName[11],"pi+ K+");
+ strcpy( fLLName[12],"pi+ proton");
+ strcpy( fLLName[13],"pi- proton");
+ strcpy( fLLName[14],"K+ K-");
+ strcpy( fLLName[15],"K+ K+");
+ strcpy( fLLName[16],"K+ proton");
+ strcpy( fLLName[17],"K- proton");
+ strcpy( fLLName[18],"deuteron deuteron");
+ strcpy( fLLName[19],"deuton alpha");
+ strcpy( fLLName[20],"triton triton");
+ strcpy( fLLName[21],"triton alpha");
+ strcpy( fLLName[22],"K0 K0");
+ strcpy( fLLName[23],"K0 K0b");
+ strcpy( fLLName[24],"deuteron triton");
+ strcpy( fLLName[25],"proton triton");
+ strcpy( fLLName[26],"proton alpha");
+ strcpy( fLLName[27],"proton lambda");
+ strcpy( fLLName[28],"neutron lambda");
+ strcpy( fLLName[29],"Lambda lambda");// gael 21May02
+ strcpy( fLLName[30],"Proton Anti-proton");// gael 21May02
+ FsiInit();
+ FsiNucl();
+
+ return *this;
+}
+
+
+double AliFemtoModelWeightGeneratorLednicky::GenerateWeight(AliFemtoPair* aPair)
+{
+ // Get hidden information pointers
+ AliFemtoModelHiddenInfo *inf1 = (AliFemtoModelHiddenInfo *) aPair->Track1()->HiddenInfo();
+ AliFemtoModelHiddenInfo *inf2 = (AliFemtoModelHiddenInfo *) aPair->Track2()->HiddenInfo();
+
+ // Calculate pair variables
+ Double_t tPx = inf1->GetTrueMomentum()->x()+inf2->GetTrueMomentum()->x();
+ Double_t tPy = inf1->GetTrueMomentum()->y()+inf2->GetTrueMomentum()->y();
+ Double_t tPz = inf1->GetTrueMomentum()->z()+inf2->GetTrueMomentum()->z();
+ Double_t tM1 = inf1->GetMass();
+ Double_t tM2 = inf2->GetMass();
+ Double_t tE1 = sqrt(tM1*tM1 + inf1->GetTrueMomentum()->mag2());
+ Double_t tE2 = sqrt(tM2*tM2 + inf2->GetTrueMomentum()->mag2());
+ Double_t tE = tE1 + tE2;
+ Double_t tPt = tPx*tPx + tPy*tPy;
+ Double_t tMt = tE*tE - tPz*tPz;//mCVK;
+ Double_t tM = sqrt(tMt - tPt);
+ tMt = sqrt(tMt);
+ tPt = sqrt(tPt);
+ Double_t tBetat = tPt/tMt;
+
+ // Boost to LCMS
+ Double_t tBeta = tPz/tE;
+ Double_t tGamma = tE/tMt;
+ fKStarLong = tGamma * (inf1->GetTrueMomentum()->z() - tBeta * tE1);
+ Double_t tE1L = tGamma * (tE1 - tBeta * inf1->GetTrueMomentum()->z());
+
+ // Rotate in transverse plane
+ fKStarOut = ( inf1->GetTrueMomentum()->x()*tPx + inf1->GetTrueMomentum()->y()*tPy)/tPt;
+ fKStarSide = (-inf1->GetTrueMomentum()->x()*tPy + inf1->GetTrueMomentum()->y()*tPx)/tPt;
+
+ // Boost to pair cms
+ fKStarOut = tMt/tM * (fKStarOut - tPt/tMt * tE1L);
+
+ tBetat = tPt/tMt;
+
+ Double_t tDX = inf1->GetEmissionPoint()->x()-inf2->GetEmissionPoint()->x();
+ Double_t tDY = inf1->GetEmissionPoint()->y()-inf2->GetEmissionPoint()->y();
+ Double_t tRLong = inf1->GetEmissionPoint()->z()-inf2->GetEmissionPoint()->z();
+ Double_t tDTime = inf1->GetEmissionPoint()->t()-inf2->GetEmissionPoint()->t();
+
+ Double_t tROut = (tDX*tPx + tDY*tPy)/tPt;
+ Double_t tRSide = (-tDX*tPy + tDY*tPx)/tPt;
+
+ fRStarSide = tRSide;
+
+ fRStarLong = tGamma*(tRLong - tBeta* tDTime);
+ Double_t tDTimePairLCMS = tGamma*(tDTime - tBeta* tRLong);
+
+ tBeta = tPt/tMt;
+ tGamma = tMt/tM;
+
+ fRStarOut = tGamma*(tROut - tBeta* tDTimePairLCMS);
+ fRStar = ::sqrt(fRStarOut*fRStarOut + fRStarSide*fRStarSide +
+ fRStarLong*fRStarLong);
+ fKStar = ::sqrt(fKStarOut*fKStarOut + fKStarSide*fKStarSide + fKStarLong*fKStarLong);
+
+ if (!SetPid(inf1->GetPDGPid(),inf2->GetPDGPid())) {
+ fWeightDen=1.;
+ return 1;
+ }
+ else { // Good Pid
+ AliFemtoThreeVector* p;
+ p=(inf1->GetTrueMomentum());
+ double p1[]={p->x(),p->y(),p->z()};
+ p=(inf2->GetTrueMomentum());
+ double p2[]={p->x(),p->y(),p->z()};
+ if ((p1[0]==p2[0])&&(p1[1]==p2[1])&&(p1[2]==p2[2])) {
+ fWeightDen=0.;
+ return 0;
+ }
+ if (fSwap) {
+ fsimomentum(*p2,*p1);
+ } else {
+ fsimomentum(*p1,*p2);
+ }
+ AliFemtoLorentzVector* tPoint;
+ tPoint=(inf1->GetEmissionPoint());
+// cout << "Pid1:dans GetWeight = " << aThPair->GetPid1() << endl;
+// cout << "Pid2:dans GetWeight = " << aThPair->GetPid2() << endl;
+// cout << "LL:in GetWeight = " << mLL << endl;
+
+ double x1[]={tPoint->x(),tPoint->y(),tPoint->z(),tPoint->t()};
+ tPoint=(inf2->GetEmissionPoint());
+ double x2[]={tPoint->x(),tPoint->y(),tPoint->z(),tPoint->t()};
+ if ((x1[0]==x2[0])&&(x1[1]==x2[1])&&(x1[2]==x2[2])&&(x1[3]==x2[3])) {
+ fWeightDen=0.;
+ return 0;
+ }
+ if (fSwap) {
+ fsiposition(*x2,*x1);
+ } else {
+ fsiposition(*x1,*x2);
+ }
+ FsiSetLL();
+ ltran12();
+ fsiw(1,fWeif,fWei,fWein);
+ if (fI3c==0) return fWein;
+ fWeightDen=fWeif;
+ return fWei;
+ };
+};
+
+AliFemtoString AliFemtoModelWeightGeneratorLednicky::Report() {
+ // create report
+ ostringstream tStr;
+ tStr << "Lednicky afterburner calculation for Correlation - Report" << endl;
+ tStr << " Setting : Quantum : " << ((fIqs) ? "On" : "Off");
+ tStr << " - Coulbomb : " << ((fIch) ? "On" : "Off") ;
+ tStr << " - Strong : " << ((fIsi) ? "On" : "Off");
+ tStr << endl;
+ tStr << " 3-Body : " << ((fI3c) ? "On" : "Off") ;
+ if (fI3c) tStr << " Mass=" << fNuclMass << " - Charge= " << fNuclCharge ;
+ tStr << endl;
+ tStr << " " << fNumProcessPair[0] << " Pairs have been Processed :" << endl;
+ int i;
+ for(i=1;i<=fLLMax;i++) {
+ if (fNumProcessPair[i])
+ tStr << " " << fNumProcessPair[i] << " " << fLLName[i] << endl;
+ }
+ if (fNumbNonId)
+ tStr << " "<< fNumbNonId << " Non Identified" << endl;
+ AliFemtoString returnThis = tStr.str();
+ return returnThis;
+}
+
+void AliFemtoModelWeightGeneratorLednicky::FsiInit(){
+ // Initialize weight generation module
+ cout << "*******************AliFemtoModelWeightGeneratorLednicky check FsiInit ************" << endl;
+ cout <<"mItest dans FsiInit() = " << fItest << endl;
+ cout <<"mIch dans FsiInit() = " << fIch << endl;
+ cout <<"mIqs dans FsiInit() = " << fIqs << endl;
+ cout <<"mIsi dans FsiInit() = " << fIsi << endl;
+ cout <<"mI3c dans FsiInit() = " << fI3c << endl;
+ fsiin(fItest,fIch,fIqs,fIsi,fI3c);
+};
+
+void AliFemtoModelWeightGeneratorLednicky::FsiNucl(){
+ // initialize weight generation taking into account the residual charge
+ cout << "*******************AliFemtoModelWeightGeneratorLednicky check FsiNucl ************" << endl;
+ cout <<"fNuclMass dans FsiNucl() = " << fNuclMass << endl;
+ cout <<"fNuclCharge dans FsiNucl() = " << fNuclCharge << endl;
+ cout <<"fNuclChargeSign dans FsiNucl() = " << fNuclChargeSign << endl;
+ fsinucl(fNuclMass,fNuclCharge*fNuclChargeSign);
+};
+
+void AliFemtoModelWeightGeneratorLednicky::FsiSetLL(){
+ // set internal pair type for the module
+ int tNS;
+ if (fSphereApp||(fLL>5)) {
+ if (fT0App) { tNS=4;}
+ else {tNS=2;};
+ } else { tNS=1;};
+ //cout <<"fLL dans FsiSetLL() = "<< fLL << endl;
+ //cout <<"tNS dans FsiSetLL() = "<< tNS << endl;
+ //cout <<"fItest dans FsiSetLL() = "<< fItest << endl;
+ llini(fLL,tNS,fItest);
+ //cout<<" end of FsiSetLL"<<endl;
+}
+
+bool AliFemtoModelWeightGeneratorLednicky::SetPid(const int aPid1,const int aPid2) {
+ // set calculated system for basing on particles' pids
+ static const int ksPi0Pid=111;
+ static const int ksPionPid=211;
+ static const int ksK0Pid=311;
+ static const int ksKPid=321;
+ static const int ksNeutPid=2112;
+ static const int ksProtPid=2212;
+ static const int ksLamPid=3122;
+ // static const int sLamLamPid=3122;
+
+ // cout << "Setting PID to " << aPid1 << " " << aPid2 << endl;
+
+ int tPidl,tPidh;
+ int tChargeFactor=1;
+
+ if (abs(aPid1)<abs(aPid2)) {
+ if (aPid1<0) tChargeFactor=-1;
+ tPidl=aPid1*tChargeFactor;
+ tPidh=aPid2*tChargeFactor;
+ fSwap=false;
+ } else {
+ if (aPid2<0) tChargeFactor=-1;
+ tPidl=aPid2*tChargeFactor;
+ tPidh=aPid1*tChargeFactor;
+ fSwap=true;
+ }
+ switch (tPidl) {
+ case ksPionPid:
+ switch (tPidh) {
+ case -ksPionPid: fLL=5; tChargeFactor*=1 ;break;
+ case ksPionPid: fLL=7; tChargeFactor*=1 ;break;
+ case -ksKPid: fLL=10;tChargeFactor*=1 ;break;
+ case ksKPid: fLL=11;tChargeFactor*=1 ;break;
+ case ksProtPid: fLL=12;tChargeFactor*=1 ;break;
+ case -ksProtPid: fLL=13;tChargeFactor*=-1;break;
+ default: fLL=0;
+ }
+ break;
+ case ksProtPid:
+ switch (tPidh) {
+ case ksProtPid: fLL=2; tChargeFactor*=1 ;break;
+ case ksLamPid: fLL=27;tChargeFactor*=1 ;break;
+ case -ksProtPid: fLL=30;tChargeFactor*=1 ;break;
+ default: fLL=0;
+ };
+ break;
+ case ksKPid:
+ switch (tPidh) {
+ case -ksKPid: fLL=14;tChargeFactor*=1 ;break;
+ case ksKPid: fLL=15;tChargeFactor*=1 ;break;
+ case ksProtPid: fLL=16;tChargeFactor*=1 ;break;
+ case -ksProtPid: fLL=17;tChargeFactor*=-1 ;break;
+ default: fLL=0;
+ };
+ break;
+ case ksK0Pid:
+ switch (tPidh) {
+ case ksK0Pid: fLL=22;tChargeFactor*=1 ;break;
+ case -ksK0Pid: fLL=23;tChargeFactor*=1 ;break;
+ default: fLL=0;
+ };
+ break;
+ case ksPi0Pid:
+ switch (tPidh) {
+ case ksPi0Pid: fLL=6; tChargeFactor*=1 ;break;
+ default: fLL=0;
+ };
+ break;
+ case ksNeutPid:
+ switch (tPidh) {
+ case ksNeutPid: fLL=1; tChargeFactor*=1 ;break;
+ case ksProtPid: fLL=3; tChargeFactor*=1 ;break;
+ case ksLamPid: fLL=28;tChargeFactor*=1 ;break;
+ default: fLL=0;
+ };
+ break; //Gael 21May02
+ case ksLamPid: //Gael 21May02
+ switch (tPidh) { //Gael 21May02
+ case ksLamPid: fLL=29;tChargeFactor*=1 ;break;//Gael 21May02
+ default: fLL=0; //Gael 21May02
+ }; //Gael 21May02
+ break; //Gael 21May02
+ default: fLL=0;
+ }
+ if (tChargeFactor!=fNuclChargeSign) {
+ fNuclChargeSign=tChargeFactor;
+ FsiNucl();
+ }
+ (fNumProcessPair[0])++;
+ if (fLL) {
+ (fNumProcessPair[fLL])++;
+ return true;
+ } else {
+ fNumbNonId++;
+ return false;
+ }
+ cout << "*******************AliFemtoModelWeightGeneratorLednicky check SetPid ************" << endl;
+ cout << "fLL=="<< fLL << endl;
+ cout << "fNuclCharge=="<< fNuclCharge << endl;
+
+}
+AliFemtoModelWeightGeneratorLednicky::~AliFemtoModelWeightGeneratorLednicky()
+{ /* no-op */ };
+
+//_____________________________________________
+void AliFemtoModelWeightGeneratorLednicky::SetPairType(Int_t aPairType)
+{
+ // set calculated system basing on the pair type
+ fPairType = aPairType;
+ if (fPairType == fgkPionPlusPionPlus) SetPid(211,211);
+ if (fPairType == fgkPionPlusPionMinus ) SetPid(211, -211);
+ if (fPairType == fgkKaonPlusKaonPlus ) SetPid(321, 321);
+ if (fPairType == fgkKaonPlusKaonMinus ) SetPid(321, -321);
+ if (fPairType == fgkProtonProton ) SetPid(2212, 2212);
+ if (fPairType == fgkProtonAntiproton ) SetPid(2212, -2212);
+ if (fPairType == fgkPionPlusKaonPlus ) SetPid(211, 321);
+ if (fPairType == fgkPionPlusKaonMinus ) SetPid(211, -321);
+ if (fPairType == fgkPionPlusProton ) SetPid(211, 2212);
+ if (fPairType == fgkPionPlusAntiproton ) SetPid(211, -2212);
+ if (fPairType == fgkKaonPlusProton ) SetPid(321, 2212);
+ if (fPairType == fgkKaonPlusAntiproton ) SetPid(321, -2212);
+}
+
+//_____________________________________________
+Int_t AliFemtoModelWeightGeneratorLednicky::GetPairType() const
+{
+ // return pair type
+ return fPairType;
+}
+
+//_____________________________________________
+void AliFemtoModelWeightGeneratorLednicky::SetPairTypeFromPair(AliFemtoPair *aPair)
+{
+ // set calculated system based on the hidden info in the pair
+ AliFemtoModelHiddenInfo *inf1 = ( AliFemtoModelHiddenInfo *) aPair->Track1()->HiddenInfo();
+ AliFemtoModelHiddenInfo *inf2 = ( AliFemtoModelHiddenInfo *) aPair->Track2()->HiddenInfo();
+
+ const Int_t ktPid1 = inf1->GetPDGPid();
+ const Int_t ktPid2 = inf2->GetPDGPid();
+
+ if (((ktPid1 == 211) && (ktPid2 == 211)) ||
+ ((ktPid1 == -211) && (ktPid2 == -211)))
+ fPairType = fgkPionPlusPionPlus;
+ else if (((ktPid1 == -211) && (ktPid2 == 211)) ||
+ ((ktPid1 == 211) && (ktPid2 == -211)))
+ fPairType = fgkPionPlusPionMinus;
+ else if (((ktPid1 == 321) && (ktPid2 == 321)) ||
+ ((ktPid1 == -321) && (ktPid2 == -321)))
+ fPairType = fgkKaonPlusKaonPlus;
+ else if (((ktPid1 == -321) && (ktPid2 == 321)) ||
+ ((ktPid1 == 321) && (ktPid2 == -321)))
+ fPairType = fgkKaonPlusKaonMinus;
+ else if (((ktPid1 == 2212) && (ktPid2 == 2212)) ||
+ ((ktPid1 == -2212) && (ktPid2 == -2212)))
+ fPairType = fgkProtonProton;
+ else if (((ktPid1 == -2212) && (ktPid2 == 2212)) ||
+ ((ktPid1 == 2212) && (ktPid2 == -2212)))
+ fPairType = fgkProtonAntiproton;
+ else if (((ktPid1 == 211) && (ktPid2 == 321)) ||
+ ((ktPid1 == -211) && (ktPid2 == -321)))
+ fPairType = fgkPionPlusKaonPlus;
+ else if (((ktPid1 == -211) && (ktPid2 == 321)) ||
+ ((ktPid1 == 211) && (ktPid2 == -321)))
+ fPairType = fgkPionPlusKaonMinus;
+ else if (((ktPid1 == 211) && (ktPid2 == 2212)) ||
+ ((ktPid1 == -211) && (ktPid2 == -2212)))
+ fPairType = fgkPionPlusProton;
+ else if (((ktPid1 == -211) && (ktPid2 == 2212)) ||
+ ((ktPid1 == 211) && (ktPid2 == -2212)))
+ fPairType = fgkPionPlusAntiproton;
+ else if (((ktPid1 == 321) && (ktPid2 == 2212)) ||
+ ((ktPid1 == -321) && (ktPid2 == -2212)))
+ fPairType = fgkKaonPlusProton;
+ else if (((ktPid1 == -321) && (ktPid2 == 2212)) ||
+ ((ktPid1 == 321) && (ktPid2 == -2212)))
+ fPairType = fgkKaonPlusAntiproton;
+ SetPid(ktPid1, ktPid2);
+}
+
+void AliFemtoModelWeightGeneratorLednicky::SetNuclCharge(const double aNuclCharge) {fNuclCharge=aNuclCharge;FsiNucl();};
+void AliFemtoModelWeightGeneratorLednicky::SetNuclMass(const double aNuclMass){fNuclMass=aNuclMass;FsiNucl();};
+
+void AliFemtoModelWeightGeneratorLednicky::SetSphere(){fSphereApp=true;};
+void AliFemtoModelWeightGeneratorLednicky::SetSquare(){fSphereApp=false;};
+void AliFemtoModelWeightGeneratorLednicky::SetT0ApproxOn(){ fT0App=true;};
+void AliFemtoModelWeightGeneratorLednicky::SetT0ApproxOff(){ fT0App=false;};
+void AliFemtoModelWeightGeneratorLednicky::SetDefaultCalcPar(){
+ fItest=1;fIqs=1;fIsi=1;fI3c=0;fIch=1;FsiInit();
+ fSphereApp=false;fT0App=false;};
+
+void AliFemtoModelWeightGeneratorLednicky::SetCoulOn() {fItest=1;fIch=1;FsiInit();};
+void AliFemtoModelWeightGeneratorLednicky::SetCoulOff() {fItest=1;fIch=0;FsiInit();};
+void AliFemtoModelWeightGeneratorLednicky::SetQuantumOn() {fItest=1;fIqs=1;FsiInit();};
+void AliFemtoModelWeightGeneratorLednicky::SetQuantumOff(){fItest=1;fIqs=0;FsiInit();};
+void AliFemtoModelWeightGeneratorLednicky::SetStrongOn() {fItest=1;fIsi=1;FsiInit();};
+void AliFemtoModelWeightGeneratorLednicky::SetStrongOff() {fItest=1;fIsi=0;FsiInit();};
+void AliFemtoModelWeightGeneratorLednicky::Set3BodyOn() {fItest=1;fI3c=1;FsiInit();FsiNucl();};
+void AliFemtoModelWeightGeneratorLednicky::Set3BodyOff() {fItest=1;fI3c=0;FsiInit();fWeightDen=1.;FsiNucl();};
+
+Double_t AliFemtoModelWeightGeneratorLednicky::GetKStar() const {return AliFemtoModelWeightGenerator::GetKStar();}
+Double_t AliFemtoModelWeightGeneratorLednicky::GetKStarOut() const { return AliFemtoModelWeightGenerator::GetKStarOut(); }
+Double_t AliFemtoModelWeightGeneratorLednicky::GetKStarSide() const { return AliFemtoModelWeightGenerator::GetKStarSide(); }
+Double_t AliFemtoModelWeightGeneratorLednicky::GetKStarLong() const { return AliFemtoModelWeightGenerator::GetKStarLong(); }
+Double_t AliFemtoModelWeightGeneratorLednicky::GetRStar() const { return AliFemtoModelWeightGenerator::GetRStar(); }
+Double_t AliFemtoModelWeightGeneratorLednicky::GetRStarOut() const { return AliFemtoModelWeightGenerator::GetRStarOut(); }
+Double_t AliFemtoModelWeightGeneratorLednicky::GetRStarSide() const { return AliFemtoModelWeightGenerator::GetRStarSide(); }
+Double_t AliFemtoModelWeightGeneratorLednicky::GetRStarLong() const { return AliFemtoModelWeightGenerator::GetRStarLong(); }
+
+AliFemtoModelWeightGenerator* AliFemtoModelWeightGeneratorLednicky::Clone() const {
+ AliFemtoModelWeightGenerator* tmp = new AliFemtoModelWeightGeneratorLednicky(*this);
+ return tmp;
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff