***** ****************************************************************** * program: bubblegum 2.0 * author: Norbert E. Ligterink, University of Pittsburgh, USA * nol1@pitt.edu, norbert.ligterink@physics.org * * date: Nov. 2002 * * synopsis: Unitarity and analytical 2x2 T-matrix from resonances * (Coherent Bubble Sum Approximation) * * changes: version 2.0 subtracts the mass shifts from the bare * resonance masses, so that the comparison between * the K-matrix and Breit-Wigner, and T-matrix is closer. * * authorship and use: use and change as you like, acknowlegdments * and citations are appreciated. * * function: calculates two-channel T-matrix for an arbitrary * number NR of resonances with mass MR(1..NR) and * complex couplings G(1..2,1..NR) to the two decay * channels with thresholds THRESH(1..2). * Instead of assuming a particular form factor, the * program allows for an expansion in (s-s_th)/s. * The threshold behavior PW(1..2) is fixed by: * Gamma_0(s) = (s-s_th)^(PW(i)/2)/s^([PW(i)/2]+1) * for channel i. * Hence PW(i) is 2l+2 in a nonrelativistic approach. * The expansion in ((s-s_th)/s) is given by: * P(s) = sum_(n=1)^NS CESP(n,i) ((s-s_th)/s)^n * Therefore the full Gamma(s) is * Gamma(s) = Gamma_0(s) P(S) = Im W(s) * The coefficients CESP(1..NS[expansion],1..2[channel]) * need to be given, or fitted from the data. * * note: The energy has to be taken as the total relativistic * energy E = SQRT(S) otherwise the result is not * covariant. * Only for the T-matrix the real part of the poles are * not fixed; they have selfenergy shifts. These shifts * can be considerable. * The singularities are regulated to make the numerical * calculations finite. Over a large range of values * the results do not depend on the regularization, * however, if the units of energy are far from the * order 1, say 10^+/-4, some deviations may occur. * In the next version, a dimensionless regularization * is planned. * In the case there is no data for the high-energy tail, * or the effects of selfenergy should be limited due * to a strong form factor, one can set the expansion * to one, with coefficients -1.00 for both channels. * * input: if INTERACTIVE = .TRUE. it requires keyboard input, * otherwise the file "in.dat" should contain the same * keystroked input (the structure is laid out in * "template.in.dat"; one number per line). * if MASSSHIFT = .TRUE. it substracts some average * mass shift from the T-matrix results, the numerical * values are given in "out.dat". Note that there is * no unambiguous mass shift definition in the coupled * channel case. * Change at lines 121-122 and 129-130 of the code. * * output: the program generates output files per channel, e.g. * t12.dat contains the 2->1 transition amplitude: * Energy (Re T_12) (Im T_12) * * files: "tij.dat" exact, unitary, analytical solution * files: "kij.dat" the K-matrix approximation * files: "bij.dat" the Breit-Wigner result (isobar model) * files: "pij.dat" the perturbative result * * the files gamma1.dat and gamma2.dat contain the real * and imaginary part of Gamma(E) = W(E); in a Hamiltonian * framework, this quantity is given by the square of * the transition matrix element: * |g_a|^2 W(E) = ||^2 * Hence the nontrivial information, to be extracted from * scattering data! * * file: out.dat contains the input data and the mass * shifts. It can be used as input file again. * ***** ****************************************************************** * dictionary * * TWO = 2 * SIXTYFOUR = 64 (max # resonances and # expansion) * * POINTS the number of points to be tabulated * STARTRANGE the starting value of the energy to be tabulated * ENDRANGE the end value of the energy to be tabulated * * NR the number of resonances (TO BE FITTED) * MR(1..NR) the (bare) masses of the resonances (TO BE FITTED) * DM(1..NR) the global mass shifts (as far as they exist) * THRESH(1..2) the threshold energies for the two decay channels * (FROM PARTICLE DATA REVIEW) * PW(1..2) the leading order behavior of Gamma(s) at threshold for * each channel (Gamma = (s-1)^(PW/2) GIVEN BY PARTIAL WAVE) * G(1..2,1..NR) the coupling constants between the decay channels * (1..2) and the resonances (1..NR) (TO BE FITTED) * NS the number of terms in the expansion of Gamma(s) * CESP(1..NS,1..2) the coefficients of the expansion of Gamma(s) * (TO BE FITTED) * W(energy,threshold energy,L,K) the real and imaginary part of * the function Im w = pi (s-s_th)^(L/2)/s^K * X(1..2,1..2) the perturbative (tree) diagram input to the LS eq. * T(1..2,1..2) the full solution (unitary and analytical) * K(1..2,1..2) the K-matrix approximate solution * B(1..2,1..2) the Breit-Wigner approximate solution * ***** ****************************************************************** PROGRAM bubblegum IMPLICIT NONE * LOGICAL INTERACTIVE, MASSSHIFT REAL*8 STARTRANGE, ENDRANGE, THRESH(2), CESP(64,2), MR(64), DM(64) INTEGER POINTS, TWO, NR, NS, PW(2), SIXTYFOUR COMPLEX*16 G(2,64) * PARAMETER(TWO = 2, SIXTYFOUR = 64) * * input from keyboard .TRUE., otherwise .FALSE. file: in.dat * * CHOOSE HERE!!!! * c INTERACTIVE = .TRUE. INTERACTIVE = .FALSE. * * perform global mass shifts to make comparison between T and K * better * * CHOOSE HERE!!!! * MASSSHIFT = .TRUE. C MASSSHIFT = .FALSE. * * IF(INTERACTIVE)THEN PRINT *,"Hello, welcome to Norbert Ligterink's bubblegum program" PRINT *,"Start energy:" READ *,STARTRANGE PRINT *,"End energy:" READ *,ENDRANGE PRINT *,"Number of points:" READ *,POINTS ELSE OPEN(UNIT=1, FILE='in.dat', STATUS='OLD') READ(1,*)STARTRANGE READ(1,*)ENDRANGE READ(1,*)POINTS ENDIF * IF(ENDRANGE .LT. 0.001 .OR. ENDRANGE .GT. 1000)THEN PRINT *,"Please, next time choose units in the order of 1." PRINT *,"It will reduce regularization errors." ENDIF * * assigning values (see there) * CALL PHYSICALVALUES(MR,G,PW,THRESH,CESP,TWO,NR,NS,SIXTYFOUR, c INTERACTIVE) * * calculation and output (in terms of energy) * generating the files:'tmatrix.dat', 'kmatrix.dat', 'perturb.dat', * 'formfactor.dat' * CALL OUTPUT(STARTRANGE,ENDRANGE,POINTS,MASSSHIFT, c MR,DM,G,PW,THRESH,CESP,TWO,NR,NS,SIXTYFOUR) * IF(.NOT. INTERACTIVE)THEN CLOSE(1) ENDIF * * dumping the input (file) and the mass shifts in "out.dat" * CALL DATAOUT(STARTRANGE,ENDRANGE,POINTS,MASSSHIFT, c MR,DM,G,PW,THRESH,CESP,TWO,NR,NS,SIXTYFOUR) * END * ***** ****************************************************************** SUBROUTINE c PHYSICALVALUES(MR,G,PW,THRESH,CESP,TWO,NR,NS,SIXTYFOUR, c INTERACTIVE) IMPLICIT NONE * * arguments * LOGICAL INTERACTIVE INTEGER TWO, SIXTYFOUR, NR, NS, PW(TWO) COMPLEX*16 G(TWO,SIXTYFOUR) REAL*8 THRESH(TWO), E, CESP(SIXTYFOUR,TWO), MR(SIXTYFOUR) * * internals * INTEGER I, J REAL*8 REG, IMG * IF(INTERACTIVE)THEN * PRINT *,"How many resonances do you have?" READ *,NR IF(NR .GT. 64)THEN PRINT *,"Error: too many, adjust the array limits" STOP ENDIF DO I=1, NR PRINT *,"The mass of resonance ",I,":" READ *,MR(I) ENDDO DO I=1, 2 DO J=1, NR PRINT *,"real part g(resonance=",J,", channel=",I,"):" READ *,REG PRINT *,"imag part g(resonance=",J,", channel=",I,"):" READ *,IMG G(I,J) = COMPLEX(REG,IMG) ENDDO ENDDO DO I=1, 2 PRINT *,"threshold energy of channel=",I,":" READ *,THRESH(I) PRINT *,"threshold behavior (partial wave) of channel=",I,":" READ *,PW(I) ENDDO * PRINT *,"Number of terms in the expansion of -Gamma(s)- :" READ *,NS IF(NS .GT. 64)THEN PRINT *,"Error: too many, adjust the array limits" STOP ENDIF DO I=1, 2 DO J=1, NS PRINT *,"the ((s-s_th)/s)^",J," coefficient of channel=",I,":" READ *,CESP(J,I) ENDDO ENDDO * ELSE * * the number of resonances * READ(1,*)NR IF(NR .GT. 64)THEN PRINT *,"Error NR: too many, adjust the array limits" STOP ENDIF * * the masses of the resonances M(1..NR[resonance number]) * DO I=1, NR READ(1,*)MR(I) ENDDO * * the coupling of the resonances to the two decay channels * G(1..2[decay channel],1..NR[resonance number]) * DO I=1, 2 DO J=1, NR READ(1,*)REG READ(1,*)IMG G(I,J) = COMPLEX(REG,IMG) ENDDO ENDDO * DO I=1, 2 * * the threshold value for the two decay channels * READ(1,*)THRESH(I) * * the "partial wave" meaning Gamma(s) goes like (s - s_th)^(PW/2) * at threshold * READ(1,*)PW(I) * ENDDO * * The number of terms in the polynomial P((s-s_th)/s) * for the expression: * Gamma(s) = (s - s_th)^(PW/2)/s^([PW/2]+1) P((s-s_th)/s) * READ(1,*)NS IF(NS .GT. 64)THEN PRINT *,"Error NS: too many, adjust the array limits" STOP ENDIF * * the coefficients in the polynomial P: * CESP(1..NS[expansion],1..2[decay channel]) * P(x) = 1 + CESP(1,channel) x + CESP(2,channel) x^2 + ... * DO I=1, 2 DO J=1, NS READ(1,*)CESP(J,I) ENDDO ENDDO * ENDIF * END * ***** ****************************************************************** SUBROUTINE DATAOUT(STARTRANGE,ENDRANGE,POINTS,MASSSHIFT, c MR,DM,G,PW,THRESH,CESP,TWO,NR,NS,SIXTYFOUR) IMPLICIT NONE * * arguments * LOGICAL INTERACTIVE, MASSSHIFT INTEGER TWO, SIXTYFOUR, NR, NS, PW(TWO), POINTS COMPLEX*16 G(TWO,SIXTYFOUR) REAL*8 THRESH(TWO), E, CESP(SIXTYFOUR,TWO), MR(SIXTYFOUR), c DM(SIXTYFOUR), STARTRANGE,ENDRANGE * * internals * INTEGER I, J * OPEN(UNIT=11, FILE='out.dat', STATUS='UNKNOWN') * WRITE(11,*)STARTRANGE," start range" WRITE(11,*)ENDRANGE," end range" WRITE(11,*)POINTS," number of points" WRITE(11,*)NR," number of resonances" * DO I=1, NR WRITE(11,*)MR(I)," mass of resonance ",I ENDDO * DO I=1, 2 DO J=1, NR WRITE(11,*)DBLE(G(I,J))," re coupling r=",J," c=",I WRITE(11,*)DIMAG(G(I,J))," im coupling r=",J," c=",I ENDDO ENDDO * DO I=1, 2 WRITE(11,*)THRESH(I)," threshold channel ",I WRITE(11,*)PW(I)," threshold behavior channel ",I ENDDO * WRITE(11,*)NS," power of the expansion of Gamma" * DO I=1, 2 DO J=1, NS WRITE(11,*)CESP(J,I)," coefficient of xi^",J," c=",I ENDDO ENDDO * IF(MASSSHIFT)THEN DO I=1, NR WRITE(11,*)DM(I)," mass shift of resonance ",I ENDDO ENDIF * * CLOSE(11) * END ***** ****************************************************************** SUBROUTINE OUTPUT(STARTRANGE,ENDRANGE,POINTS,MASSSHIFT, c MR,DM,G,PW,THRESH,CESP,TWO,NR,NS,SIXTYFOUR) IMPLICIT NONE * * arguments * LOGICAL MASSSHIFT REAL*8 STARTRANGE, ENDRANGE INTEGER POINTS, TWO, SIXTYFOUR, NR, NS, PW(TWO) COMPLEX*16 G(TWO,SIXTYFOUR) REAL*8 THRESH(TWO), CESP(SIXTYFOUR,TWO), MR(SIXTYFOUR), c DM(SIXTYFOUR) * * internals * INTEGER I REAL*8 E COMPLEX*16 T(TWO,TWO),K(TWO,TWO),X(TWO,TWO),B(TWO,TWO),WW(TWO) * OPEN(UNIT=21, FILE='t11.dat', STATUS='UNKNOWN') OPEN(UNIT=22, FILE='t12.dat', STATUS='UNKNOWN') OPEN(UNIT=23, FILE='t21.dat', STATUS='UNKNOWN') OPEN(UNIT=24, FILE='t22.dat', STATUS='UNKNOWN') OPEN(UNIT=31, FILE='k11.dat', STATUS='UNKNOWN') OPEN(UNIT=32, FILE='k12.dat', STATUS='UNKNOWN') OPEN(UNIT=33, FILE='k21.dat', STATUS='UNKNOWN') OPEN(UNIT=34, FILE='k22.dat', STATUS='UNKNOWN') OPEN(UNIT=41, FILE='p11.dat', STATUS='UNKNOWN') OPEN(UNIT=42, FILE='p12.dat', STATUS='UNKNOWN') OPEN(UNIT=43, FILE='p21.dat', STATUS='UNKNOWN') OPEN(UNIT=44, FILE='p22.dat', STATUS='UNKNOWN') OPEN(UNIT=71, FILE='b11.dat', STATUS='UNKNOWN') OPEN(UNIT=72, FILE='b12.dat', STATUS='UNKNOWN') OPEN(UNIT=73, FILE='b21.dat', STATUS='UNKNOWN') OPEN(UNIT=74, FILE='b22.dat', STATUS='UNKNOWN') OPEN(UNIT=51, FILE='gamma1.dat', STATUS='UNKNOWN') OPEN(UNIT=52, FILE='gamma2.dat', STATUS='UNKNOWN') * * calculating the massshifts DM(1..NR) for the T-matrix results * IF(MASSSHIFT)THEN CALL MASSHIFT(STARTRANGE,ENDRANGE,POINTS, c MR,DM,G,PW,THRESH,CESP,TWO,NR,NS,SIXTYFOUR) ENDIF * * DO I=0, POINTS-1 * E = STARTRANGE + (ENDRANGE-STARTRANGE)*DBLE(I)/DBLE(POINTS-1) * CALL TMATRIX(T,K,X,B,WW,PW,TWO,THRESH,G,MR,DM,NR,E,CESP,NS,64, c MASSSHIFT,0,0.0D0) * * full T matrix * WRITE(21,222),E,DBLE(T(1,1)),DIMAG(T(1,1)) WRITE(22,222),E,DBLE(T(1,2)),DIMAG(T(1,2)) WRITE(23,222),E,DBLE(T(2,1)),DIMAG(T(2,1)) WRITE(24,222),E,DBLE(T(2,2)),DIMAG(T(2,2)) * * the K-matrix (no real parts, i.e., no analyticity) * WRITE(31,222),E,DBLE(K(1,1)),DIMAG(K(1,1)) WRITE(32,222),E,DBLE(K(1,2)),DIMAG(K(1,2)) WRITE(33,222),E,DBLE(K(2,1)),DIMAG(K(2,1)) WRITE(34,222),E,DBLE(K(2,2)),DIMAG(K(2,2)) * * perturbative result (tree diagram) * WRITE(41,222),E,DBLE(X(1,1)),DIMAG(X(1,1)) WRITE(42,222),E,DBLE(X(1,2)),DIMAG(X(1,2)) WRITE(43,222),E,DBLE(X(2,1)),DIMAG(X(2,1)) WRITE(44,222),E,DBLE(X(2,2)),DIMAG(X(2,2)) * * breit-wigner result (isobar diagrams; tree with widths) * WRITE(71,222),E,DBLE(B(1,1)),DIMAG(B(1,1)) WRITE(72,222),E,DBLE(B(1,2)),DIMAG(B(1,2)) WRITE(73,222),E,DBLE(B(2,1)),DIMAG(B(2,1)) WRITE(74,222),E,DBLE(B(2,2)),DIMAG(B(2,2)) * * the real and imaginary part of Gamma(E) = W(E) * WRITE(51,222),E,DIMAG(WW(1)),DBLE(WW(1)) WRITE(52,222),E,DIMAG(WW(2)),DBLE(WW(2)) * ENDDO * CLOSE(21) CLOSE(22) CLOSE(23) CLOSE(24) CLOSE(31) CLOSE(32) CLOSE(33) CLOSE(34) CLOSE(41) CLOSE(42) CLOSE(43) CLOSE(44) CLOSE(71) CLOSE(72) CLOSE(73) CLOSE(74) CLOSE(51) CLOSE(52) * DO I=1, NR C PRINT *,"Resonance:",I," mass=",MR(I)," mass shift=",DM(I) ENDDO * 222 FORMAT(9G14.5) * END * ***** ****************************************************************** * * the analytic function W of the dimension of energy which * imaginary part is given by: * * 2 2 L/2 * PI THRESH (E - THRESH ) * Im W = --------------------------- * E 2k * (-------) * THRESH * * ***** ****************************************************************** COMPLEX*16 FUNCTION W(E,THRESHOLD,L,K) IMPLICIT NONE REAL*8 E, THRESHOLD, S, WR, WI, POL, COEFF, LOGCOEFF REAL*8 PI, ATAN, DLOG, DEXP, SIG INTEGER L, K, I, J * PI = 4.0D0*ATAN(1.0D0) * * conversion to a threshold normalised S variable * S = E*E/(THRESHOLD*THRESHOLD) * * avoiding floating point errors at threshold and the origin * IF(ABS(S-1.0D0) .LT. 0.000001D0)THEN S = 0.999999D0 ENDIF * IF(ABS(S) .LT. 0.0001D0)THEN S = 0.0001D0 ENDIF * * the imaginary part * IF(S .GT. 1.0D0)THEN WI = PI*DEXP(0.5D0*DBLE(L)*DLOG(S-1.0D0))/DEXP(DBLE(K)*DLOG(S)) ELSE WI = 0.0D0 ENDIF * * the real part * * square root formula for L is odd * IF(MOD(L,2) .NE. 0)THEN SIG = (-1.0D0)**((L-1)/2) IF(S .LT. 1.0D0)THEN WR = SIG*PI*DEXP(0.5D0*DBLE(L)*DLOG(1.0D0-S))/ c DEXP(DBLE(K)*DLOG(S)) ELSE WR = 0.0D0 ENDIF * * the coefficients of the series expansion are given by * generalized binomial coefficients * POL = -SIG*PI*DEXP(-DBLE(K)*DLOG(S)) DO I=1, K C PRINT *,POL," S**",-K-1+I WR = WR + POL POL = -POL*S*DBLE(L-2*I+2)/DBLE(2*I) ENDDO * * log formula for L is even * ELSE SIG = (-1.0D0)**(L/2) IF(S .LT. 1.0D0)THEN WR = SIG*DLOG(1.0D0-S)*DEXP(0.5D0*DBLE(L)* c DLOG(1.0D0-S))/DEXP(DBLE(K)*DLOG(S)) ELSE WR = DLOG(S-1.0D0)*DEXP(0.5D0*DBLE(L)* c DLOG(S-1.0D0))/DEXP(DBLE(K)*DLOG(S)) ENDIF * * starting at 1/s^k, the Taylor expansion is generated * rescursively, term by term COEFF is the coefficient * of (1-S)^power, LOGCOEFF is the coefficient of * (1-S)^power LOG|1-S| * POL = -SIG*DEXP(-DBLE(K)*DLOG(S)) COEFF = 0.0D0 LOGCOEFF = 1.0D0 DO I=1, K C PRINT *,-COEFF," S**",-K-1+I WR = WR + POL*COEFF COEFF = -(0.5D0*COEFF*DBLE(L-2*I+2) + LOGCOEFF)/DBLE(I) LOGCOEFF = -0.5D0*LOGCOEFF*DBLE(L-2*I+2)/DBLE(I) POL = S*POL ENDDO ENDIF * * combining real and imaginary part (COMPLEX*16) * (COMPLEX concatenates only, at least for GNU fortran) * W = COMPLEX(THRESHOLD*WR,THRESHOLD*WI) * END * ***** ****************************************************************** SUBROUTINE TMATRIX(T,K,X,B,WW,PW,TWO,THRESH,G, c MR,DM,NR,E,CESP,NS,SIXTYFOUR,MASSSHIFT,RES,DELTA) IMPLICIT NONE * * arguments * LOGICAL MASSSHIFT INTEGER TWO, SIXTYFOUR, NR, NS, PW(TWO), RES COMPLEX*16 T(TWO,TWO), K(TWO,TWO), X(TWO,TWO), B(TWO,TWO) COMPLEX*16 G(TWO,SIXTYFOUR), WW(TWO) REAL*8 THRESH(TWO), E, CESP(SIXTYFOUR,TWO), MR(SIXTYFOUR), c DM(SIXTYFOUR), DELTA * * internals * INTEGER I, J, L COMPLEX*16 FAC, DETX, DENOMINATOR, W, WK(2) REAL*8 EPS, DDM(64) * * the regularization of the perturbative singularity. * EPS = 0.000001D0 * IF(MASSSHIFT)THEN DO I=1, NR DDM(I) = DM(I) ENDDO ELSE DO I=1, NR DDM(I) = 0.0D0 ENDDO DDM(RES) = DELTA ENDIF * * * initializing X: the perturbative input * DO I=1, 2 DO J=1, 2 X(I,J) = COMPLEX(0.0D0,0.0D0) ENDDO ENDDO * * if NS=0 only the leading term participates, otherwise the * coefficients CESP(1..NS) need to be set * DO J=1, 2 WW(J) = W(E,THRESH(J),PW(J),PW(J)/2+1) DO I=1, NS WW(J) = WW(J) + COMPLEX(CESP(I,J),0.0D0)* c W(E,THRESH(J),PW(J)+2*I,PW(J)/2+I+1) ENDDO IF(DIMAG(WW(J)) .LT. -0.0000001D0)THEN PRINT *,"The expansion coefficients are too large:" PRINT *,"Gamma(channel=",J,") is negative" PRINT *,"exciting" STOP ENDIF ENDDO * * * a small complex part removes the singularity, but is * of no further consequence since the singularity * cancels between the numerator and the denominator of T * The maximum violations of unitarity are typically of the * order of ten times this number: * |Im(T_ij) - sum_k T^*_ik T_kj| < 10*0.000001D0 * * the mass shifts DM are subtracted for the T-matrix results * making the bare resonance mass the physical mass, approximately * DO L=1, NR FAC = COMPLEX(1.0D0,0.0D0)/(COMPLEX(E,0.0D0) - c COMPLEX(MR(L)-DDM(L),EPS)) DO I=1, 2 DO J=1, 2 X(I,J) = X(I,J) + G(I,L)*CONJG(G(J,L))*FAC ENDDO ENDDO ENDDO * DETX = X(1,1)*X(2,2) - X(1,2)*X(2,1) * * the common denominator for all the T * DENOMINATOR = COMPLEX(1.0D0,0.0D0) - c X(1,1)*WW(1) - X(2,2)*WW(2) + DETX*WW(1)*WW(2) * IF(.NOT. MASSSHIFT)THEN DM(64) = DBLE((DENOMINATOR-1.0D0)*(COMPLEX(E,0.0D0) c -COMPLEX(MR(RES)-DDM(RES),EPS))) + E - MR(RES) + DDM(RES) ENDIF * * reducing the superfluous calculations for the masshift routine * IF(RES .NE. 0)THEN GOTO 999 ENDIF * * the reduced transition amplitudes * T(1,1) = (X(1,1) - DETX*WW(2))/DENOMINATOR * T(1,2) = X(1,2)/DENOMINATOR * T(2,1) = X(2,1)/DENOMINATOR * T(2,2) = (X(2,2) - DETX*WW(1))/DENOMINATOR * * setting the real part and the mass shifts to zero for * an analogous K-matrix result * DO I=1, 2 DO J=1, 2 X(I,J) = COMPLEX(0.0D0,0.0D0) ENDDO ENDDO * DO L=1, NR FAC = COMPLEX(1.0D0,0.0D0)/(COMPLEX(E,0.0D0) - c COMPLEX(MR(L),0.000001D0)) DO I=1, 2 DO J=1, 2 X(I,J) = X(I,J) + G(I,L)*CONJG(G(J,L))*FAC ENDDO ENDDO ENDDO * DETX = X(1,1)*X(2,2) - X(1,2)*X(2,1) * DO I=1, 2 WK(I) = COMPLEX(0.0D0,DIMAG(WW(I))) ENDDO * * the common denominator for all the K * DENOMINATOR = COMPLEX(1.0D0,0.0D0) - c X(1,1)*WK(1) - X(2,2)*WK(2) + DETX*WK(1)*WK(2) * * the reduced transition amplitudes * K(1,1) = (X(1,1) - DETX*WK(2))/DENOMINATOR * K(1,2) = X(1,2)/DENOMINATOR * K(2,1) = X(2,1)/DENOMINATOR * K(2,2) = (X(2,2) - DETX*WK(1))/DENOMINATOR * * the Breit-Wigner approximation (the isobar model) * DO I=1, 2 DO J=1, 2 B(J,I) = COMPLEX(0.0D0,0.0D0) DO L=1, NR B(J,I) = B(J,I) + c CONJG(G(J,L))*G(I,L)/(E - MR(L) - c CONJG(G(1,L))*G(1,L)*WK(1) - CONJG(G(2,L))*G(2,L)*WK(2)) ENDDO ENDDO ENDDO * * the full transition amplitudes ("reduced T" times "phase space") * (the imaginary parts are the same for T, X, B and K) * Flipping the sign on the real part: conventions, conventions. * DO I=1, 2 DO J=1, 2 FAC = SQRT(DIMAG(WW(I))*DIMAG(WW(J))) X(I,J) = -CONJG(X(I,J))*FAC T(I,J) = -CONJG(T(I,J))*FAC K(I,J) = -CONJG(K(I,J))*FAC B(I,J) = -CONJG(B(I,J))*FAC ENDDO ENDDO * * from the mass shift calculation * 999 CONTINUE * END ***** ****************************************************************** * the mass shifts are determined numerically by taking the common * denominator close to E = MR(I). In that case X is dominated * by the contribution of the I-th resonance. The left shift DM(I)<0 * and the right shift DM(I)>0 are combined to reduce the numerical * dependence on the regularization of the singularity in X. * The real part of the denominator DM(I) in an expansion around * (E-MR(I)) divided by the derivative and is subtracted from MR(I). * SUBROUTINE MASSHIFT(STARTRANGE,ENDRANGE,POINTS, c MR,DM,G,PW,THRESH,CESP,TWO,NR,NS,SIXTYFOUR) IMPLICIT NONE * * arguments * REAL*8 STARTRANGE, ENDRANGE INTEGER POINTS, TWO, SIXTYFOUR, NR, NS, PW(TWO) COMPLEX*16 G(TWO,SIXTYFOUR) REAL*8 THRESH(TWO), CESP(SIXTYFOUR,TWO), MR(SIXTYFOUR), c DM(SIXTYFOUR) * * internals * INTEGER I REAL*8 E, RIGHT, LEFT, DELTA COMPLEX*16 T(TWO,TWO),K(TWO,TWO),X(TWO,TWO),B(TWO,TWO),WW(TWO) * * DM(I) = +/-DELTA, the shift should be smaller than resolution * of the plot, but not too small that it suffers from errors due * to the regularization of X, i.e., the perturbative input with * zero (= 0.00001D0) width. * DELTA = 0.001D0*MR(1) * DO I=1, NR * E = MR(I) * * Better to have one-and-the-same routine doing all calculations * CALL TMATRIX(T,K,X,B,WW,PW,TWO,THRESH,G,MR,DM,NR,E,CESP,NS,64, c .FALSE.,I,DELTA) * * DM(64) contains the common denominator * LEFT = DM(64) * * CALL TMATRIX(T,K,X,B,WW,PW,TWO,THRESH,G,MR,DM,NR,E,CESP,NS,64, c .FALSE.,I,-DELTA) * RIGHT = DM(64) * * in the odd combination the imaginary part from the * regularization drops out in leading order. * * rather arbitrary set of rules to exclude unrealistic massshifts * IF(ABS(1.0D0 + 0.5D0*(RIGHT - LEFT)/DELTA) .LT. 0.9D0 .AND. c ABS(DELTA*(RIGHT+LEFT)/(RIGHT - LEFT)) .LT. 0.80D0*MR(I))THEN DM(I) = DELTA*(RIGHT+LEFT)/(RIGHT-LEFT) ELSE PRINT *,"failed to find a appropriate mass shift for r=",I DM(I) = 0.0D0 ENDIF * ENDDO * * END ***** ****************************************************************** * end bubblegum ***** ******************************************************************