SetDisplayInfo added
[u/mrichter/AliRoot.git] / HBTAN / fsiini.F
1
2           SUBROUTINE FSIINI
3 C---Note:
4 C-- ICH= 0 (1) if the Coulomb interaction is absent (present);
5 C-- ISPIN= JJ= 1,2,..,MSPIN denote increasing values of the pair
6 C-- total spin S.
7 C-- To calculate the CF of two particles (with masses m1, m2 and
8 C-- charges C1, C2) the following information is required:
9 C-- AM= twice the reduced mass= 2*m1*m2/(m1+m2) in GeV/c^2,
10 C-- DM= (m1-m2)/(m1+m2), required if NS=2;
11 C-- AC= Bohr radius= 2*137.036*0.1973/(C1*C2*AMH) in fm;
12 C-- AC > 1.D9 if C1*C2= 0, AC < 0 if C1*C2 < 0;
13 C-- MSPIN= MSPINH(LL)= number of the values of the total pair spin S;
14 C-- FD= FDH(LL,JJ), RD= RDH(LL,JJ)= scattering length and effective
15 C-- radius for each value of the total pair spin S, JJ= 1,..,MSPIN;     ;
16 C-- the corresponding square well parameters EB= EBH(LL,JJ), RB=
17 C-- RBH(LL,JJ) (required if NS=1) may be calculated by SEAR;
18 C-- if the effective range approximation is not valid (as is the case,
19 C-- e.g., for two-pion system) a code for calculation of the scattering
20 C-- amplitude should be supplemented;
21 C-- RHO= RHOH(LL,JJ), SF= SFH(LL,JJ), SE= SEH(LL) are spin factors;
22 C-- RHO= the probability that the spins j1 and j2 of the two particles
23 C-- will combine in a total spin S;
24 C-- RHO= (2*S+1)/[(2j1+1)*(2j2+1)] for unpolarized particles;
25 C-- RHO= (1-P1*P2)/4 and (3+P1*P2)/4 correspond to S=0 and 1 in the
26 C-- case of spin-1/2 particles with polarizations P1 and P2;
27 C-----------------------------------------------------------------------
28       IMPLICIT REAL*8 (A-H,O-Z)
29       COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
30       COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS,
31      1               X,Y,Z,T,RP,RPS
32       COMMON/FSI_SPIN/RHO(10)
33       COMMON/FSI_ACH/HPR,AC,ACH,ACHR,ISPIN,MSPIN
34       COMMON/FSI_NS/LL,NS,ICH,ISI,IQS,I3C,I3S
35       COMMON/FSI_FD/FD(10),RD(10)
36       COMMON/FSI_C/C(10),AM,AMS,DM
37       COMMON/FSI_CONS/PI,PI2,SPI,DR,W
38       COMPLEX*16 C
39       COMMON/FSI_AA/AA
40       COMMON/FSI_AAPI/AAPI(20,2)/FSI_AAND/AAND(20,4)
41       COMMON/FSI_SW/RB(10),EB(10),BK(10),CDK(10),SDK(10),
42      1              SBKRB(10),SDKK(10)
43
44       COMMON/LEDWEIGHT/WEIF,WEI,WEIN,ITEST,IRANPOS 
45
46       DIMENSION FDH(30,10),RDH(30,10),EBH(30,10),RBH(30,10)
47       DIMENSION RHOH(30,10)
48       DIMENSION AM1H(30),AM2H(30),C1H(30),C2H(30),MSPINH(30)
49 C============= declarations pour l'appel de READ_FILE()============
50       CHARACTER*10 KEY
51       CHARACTER*8  CH8
52       INTEGER*4    INT4
53       REAL*8       REAL8
54       INTEGER*4    IERR
55 C
56 C--- mass of the first and second particle
57       DATA AM1H/.9395656D0,.9382723D0,.9395656D0,3.7294D0,.13957D0,
58      C          .13498D0,.13957D0, .9395656D0, .9382723D0,
59      C          4*.13957D0,4*.493677D0,
60      C          2*1.875613D0,2*2.808D0,2*.497672D0,
61      C          1.875613D0,2*.9382723D0, 4*0.D0/
62       DATA AM2H/.9395656D0,.9382723D0,.9382723D0,3.7294D0,.13957D0,
63      C          .13498D0,.13957D0, 2*1.875613D0,
64      C          2*.493677D0,2*.9382723D0,
65      C          2*.493677D0,2*.9382723D0,
66      C          1.875613D0,3.7294D0,2.808D0,3.7294D0,
67      C          2*.497672D0,2*2.808D0,3.7294D0, 4*0.D0/
68 C---  charge of the first and second particle
69       DATA C1H/0.D0,1.D0,0.D0,2.D0, 1.D0,0.D0,1.D0,0.D0,1.D0,
70      C         3*1.D0,-1.D0,3*1.D0,-1.D0,
71      C         4*1.D0,2*0.D0,3*1.D0, 4*0.D0/
72       DATA C2H/0.D0,1.D0,1.D0,2.D0,-1.D0,0.D0,3*1.D0,
73      C         -1.D0,3*1.D0,-1.D0,3*1.D0,
74      C         1.D0,2.D0,1.D0,2.D0,2*0.D0,2*1.D0,2.D0, 4*0.D0/
75 C---MSPIN vs (LL)
76       DATA MSPINH/3*2,4*1,2*2,8*1,3,1,2,1,2*1,2*2,1, 4*0/
77 C---Spin factors RHO vs (LL,ISPIN)
78       DATA RHOH/3*.25D0, 4*1.D0, 2*.3333D0, 8*1.D0, 
79      C          .1111D0,1.D0,.25D0,1.D0,2*1.D0,
80      1          .3333D0,.25D0,1.D0, 4*0.D0,
81      C          3*.75D0, 4*0.D0, 2*.6667D0, 8*0.D0, 
82      C          .3333D0,.0D0,.75D0,.0D0,2*0.D0,
83      2          .6667D0,.75D0,0.D0, 4*0.D0, 
84      C          17*.0D0,.5556D0,3*0.D0, 5*0.D0,4*0.D0,210*0.D0/
85 C---Scattering length FD and effective radius RD in fm vs (LL,ISPIN)
86       DATA FDH/17.0D0,7.8D0,23.7D0,2230.1218D0,.225D0,.081D0,-.063D0,
87      C        -.65D0,-2.73D0,
88      C        .137D0,-.071D0,-.148D0,.112D0,2*1.D-6,-.360D0,
89      1        2*1.D-6,1.344D0,6*1.D-6,-5.628D0, 4*0.D0,
90      C -10.8D0,2*-5.4D0,4*0.D0,-6.35D0,-11.88D0,8*0.D0,9*0.D0,4*0.D0,
91      C 240*0.D0/
92       DATA RDH/2.7D0,2.8D0,2.7D0,1.12139906D0,-44.36D0,64.0D0,784.9D0,
93 c--------|---------|---------|---------|---------|---------|---------|----------      
94      C  477.9D0, 2.27D0, 9*0.D0,-69.973D0, 6*0.D0,3.529D0, 4*0.D0,
95      C       3*1.7D0,4*0.D0,2.0D0,2.63D0, 17*0.D0, 4*0.D0, 240*0.D0/
96 C---Corresponding square well parameters RB (width in fm) and
97 C-- EB =SQRT(-AM*U) (in GeV/c); U is the well height
98       DATA RBH/2.545739D0,   2.779789D0, 2.585795D0, 5.023544D0,
99      C .124673D0, .3925180D0,.09D0, 2.D0, 4.058058D0, 17*0.D0, 4*0.D0,
100      C         3*2.003144D0,
101      C         4*0.D0, 2.D0, 4.132163D0, 17*0.D0, 4*0.D0, 240*0.D0/
102       DATA EBH/.1149517D0,    .1046257D0,   .1148757D0, .1186010D0,
103      C    .7947389D0,2.281208D0,8.7D0,.4D0,.1561219D0,17*0.D0,4*0.D0,
104      C         3*.1847221D0,
105      C         4*0.D0, .4D0, .1150687D0, 17*0.D0, 4*0.D0, 240*0.D0/
106 C=======< constants >========================
107       W=1/.1973D0    ! from fm to 1/GeV
108       PI=4*DATAN(1.D0)
109       PI2=2*PI
110       SPI=DSQRT(PI)
111       DR=180.D0/PI   ! from radian to degree
112         AC1=1.D10
113         AC2=1.D10
114  
115 C---setting particle masses and charges
116       AM1=AM1H(LL)
117       AM2=AM2H(LL)
118       C1=C1H(LL)
119       C2=C2H(LL)
120  
121 C-- Switches:
122 C   ISI=1(0)  the strong interaction between the two particles ON (OFF)
123 C   IQS=1(0)  the quantum statistics ON (OFF);
124 C             should be OFF for nonidentical particles
125 C   I3C=1(0)  the Coulomb interaction with the nucleus ON (OFF)
126 C   I3S=1(0)  the strong interaction with the nucleus ON (OFF)
127 C   ICH=1(0)  if C1*C2 is different from 0 (is equal to 0)
128 C-  to switch off the Coulomb force between the two particles
129 C   put ICH=0 and substitute the strong amplitude parameters by
130 C   the ones not affected by Coulomb interaction
131  
132       IF(ITEST.EQ.0)THEN
133
134        ICH=0
135        IF(C1*C2.NE.0.D0) ICH=1
136        IQS=0
137        IF(C1+AM1.EQ.C2+AM2) IQS=1
138        I3S=0  ! only this option is available
139        ISI=1
140        I3C=1
141        
142        ENDIF
143        
144 c23456       
145        write(*,*)'FSIINI ITEST ich iqs i3s isi i3c',ITEST, 
146      +  ICH,IQS,I3S,ISI,I3C
147  
148 C==================================================================
149 C---fm to 1/GeV
150       DO 3 J1=1,30
151       DO 3 J2=1,10
152       FDH(J1,J2)=FDH(J1,J2)*W
153       RDH(J1,J2)=RDH(J1,J2)*W
154  3    RBH(J1,J2)=RBH(J1,J2)*W
155 C---calcul. twice the reduced mass (AM), the relative mass difference
156 C-- (DM) and the Bohr radius (AC)
157       AM=2*AM1*AM2/(AM1+AM2)
158       AMS=AM*AM
159       DM=(AM1-AM2)/(AM1+AM2)
160       AC=1.D10
161       C12=C1*C2
162       IF(C12.NE.0.D0)AC=2*137.036D0/(C12*AM)
163 C---Setting spin factors
164       MSPIN=MSPINH(LL)
165       MSP=MSPIN
166       DO 91 ISPIN=1,10
167   91  RHO(ISPIN)=RHOH(LL,ISPIN)
168 C---Integration limit AA in the spherical wave approximation
169       AA=0.D0
170 cc      IF(NS.EQ.2.OR.NS.EQ.4)AA=.5D0 !!in 1/GeV --> 0.1 fm
171       IF(NS.EQ.2.OR.NS.EQ.4)AA=6.D0 !!in 1/GeV --> 1.2 fm
172 C---Setting scatt. length (FD), eff. radius (RD) and, if possible,
173 C-- also the corresp. square well parameters (EB, RB)
174       DO 55 JJ=1,MSP
175       ISPIN=JJ
176       FD(JJ)=FDH(LL,JJ)
177       RD(JJ)=RDH(LL,JJ)
178       EB(JJ)=EBH(LL,JJ)
179       RB(JJ)=RBH(LL,JJ)
180       IF(LL.NE.8.AND.LL.NE.9)GOTO25
181 C---Resets FD and RD for a nucleon-deuteron system (LL=8,9)
182       JH=LL-7+2*JJ-2
183       FD(JJ)=AAND(1,JH)
184       RD(JJ)=AAND(2,JH)-2*AAND(3,JH)/AAND(1,JH)
185 C---Resets FD and RD for a pion-pion system (LL=5,6,7)
186   25  IF(LL.GT.7.OR.LL.LT.5)GOTO 24
187       IF(LL.EQ.7)FD(JJ)=AAPI(1,2)/AM
188       IF(LL.EQ.5)FD(JJ)=(.6667D0*AAPI(1,1)+.3333D0*AAPI(1,2))/AM
189       IF(LL.EQ.6)FD(JJ)=(.3333D0*AAPI(1,1)+.6667D0*AAPI(1,2))/AM
190       AKS=0.D0
191       DAKS=1.D-5
192       AKSH=AKS+DAKS
193       AKH=DSQRT(AKSH)
194       GPI1H=GPIPI(AKSH,1)
195       GPI2H=GPIPI(AKSH,2)
196       H=1/FD(JJ)
197       IF(LL.EQ.7)C(JJ)=1/DCMPLX(GPI2H,-AKH)
198       IF(LL.EQ.5)
199      +C(JJ)=.6667D0/DCMPLX(GPI1H,-AKH)+.3333D0/DCMPLX(GPI2H,-AKH)
200       IF(LL.EQ.6)
201      +C(JJ)=.3333D0/DCMPLX(GPI1H,-AKH)+.6667D0/DCMPLX(GPI2H,-AKH)
202       HH=DREAL(1/C(JJ))
203       RD(JJ)=2*(HH-H)/DAKS
204  24   CONTINUE
205 C---Calculation continues for any system (any LL)
206  55   CONTINUE
207       END
208 C=======================================================
209
210
211
212
213
214
215
216
217
218
219
220       
221        FUNCTION GPIPIold(X,J)                                                                    
222 C--- GPIPI = k*COTG(DELTA), X=k^2                                                            
223 C--  J=1(2) corresponds to isospin=0(2)                                                      
224        IMPLICIT REAL*8 (A-H,O-Z)                                                              
225 c--       Include 'common_fsi_aapi.inc'                                                          
226 c--       Include 'common_fsi_c.inc'                                                             
227        COMMON/FSI_AAPI/AAPI(20,2)                                                         
228        COMMON/FSI_C/HELP(20),AM,AMS,DM                                                    
229        OM=DSQRT(X+AMS)                                                                        
230        XX=X/AMS                                                                               
231        GPIPI=OM/AAPI(1,J)                                                                     
232        GPIPI=GPIPI*(1+(AAPI(3,J)-AAPI(1,J)**2)*XX+AAPI(4,J)*XX*XX)                            
233        GPIPI=GPIPI/(1+(AAPI(3,J)+AAPI(2,J)/AAPI(1,J))*XX)                                     
234        GPIPIOLD=GPIPI                                                                
235        END   
236       
237        FUNCTION GPIN(X,J)                                                                    
238 C--- GPIN = k*COTG(DELTA), X=k^2                                                            
239 C--  J=1(2) corresponds to piN isospin=1/2(3/2)                                             
240        IMPLICIT REAL*8 (A-H,O-Z)                                                             
241 c--       Include 'common_fsi_aapin.inc'                                                        
242        COMMON/FSI_AAPIN/AAPIN(20,2)                                                      
243        GPIN=1/AAPIN(1,J)+.5D0*AAPIN(2,J)*X                                                   
244                                                                                        
245        END