DEMPgen/include/PhaseShift.h

0001 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0002   This file is the C++ adaptation of a Fortran routine for SIMC,
0003   parameterizing the cross-sections for pion-nucleon elastic scattering.
0004   The adaptation was done by Zafar Ahmed for use in the original DEMP event
0005   generator and is included here unmodified.
0006   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
0007
0008
0009 #include <stdio.h>
0010 #include <string>
0011 #include <iostream>
0012 #include <fstream>
0013 #include <vector>
0014 #include <algorithm>
0015 #include <map>
0016 #include <cmath>
0017 #include <iomanip>
0018 #include <fstream>
0019 #include <sstream>
0020 #include <cstdlib>
0021 #include <complex>
0022
0023 using namespace std;
0024 using std::setw;
0025 using std::setprecision;
0026 using std::cout;
0027 using std::cin;
0028 using std::endl;
0029
0030 double fZ0 = 0;
0031 double fZ1 = 0;
0032 double fZ2 = 0;
0033
0034
0035 /* ************************************************** */
0036 /* ! written by A. Shinozaki, Oct/00. */
0037
0038 /* ! this function subroutine returns equation 2.25 (p. 812) of Feshbach, */
0039 /* ! "Theoretical Nuclear Physics".  This is the parameterization of the */
0040 /* ! pi-N phase shifts for pion energies less than 400 MeV (LAB). */
0041
0042 /* ! s = pi-N CM energy^2 */
0043 /* ! k = pi-N CM momentum */
0044 /* ! l = pi-N CM angular momentum */
0045 /* ! ij= parameter selecting which set of phaseshift parameters to use */
0046 double formula( double s , double k, int l, int ij ) {
0047
0048   ij = ij - 1;
0049
0050   /*   ********************************************************************** */
0051   /* * Phase shift data are from G. Rowe, M. Solomon and R.H. Landau */
0052   /* ! PRC 18 (1978) 584. */
0053   /* ! L_2t2j=          S11      S31      P11       P13      P31      P33 */
0054   double x[6]   = {    0.44,    0.31,    0.61,      0.23,    0.22,    0.99 };
0055   double ss0[6] = {   1.550,   1.655,   1.435,     1.815,   1.850,   1.233 };
0056   double k0[6]  = {   0.477,   0.550,   0.393,     0.656,   0.678,   0.228 };
0057   double G0[6]  = {   0.105,   0.170,   0.230,     0.255,   0.200,   0.116 };
0058   double b[6]   = {   0.168,  -0.112,  -0.0571,  -0.0131, -0.0291,   0.114 };
0059   double c[6]   = { -0.0354, -0.0307,   0.0254,  0.00122, 0.00345, -0.0154 };
0060   double d[6]   = {  27.E-4,  21.E-4,  -29.E-4,  -0.4E-4, -1.5E-4,  7.2E-4 };
0061
0062   double xmpi = 0.13957;
0063   double kom = k/xmpi;
0064   double ss = s;
0065   double fFormula = 0;
0066
0067   if ( k > 0.353 ){
0068     kom = 0.353 / xmpi;
0069     ss  = pow( 1.383 , 2);
0070   }
0071
0072   double kom2 = pow( kom , 2);
0073   int l21  = 2 * l + 1;
0074   double s0 = pow( ss0[ ij ] , 2 );
0075
0076   if ( std::abs( s0 - ss ) < 1.e-10 ) {
0077     fFormula = 1.0;
0078   }
0079   else {
0080     double tangent = ( ( pow( kom , l21 ) ) *
0081                ( b[ ij ] + ( c[ ij ] + d[ ij ] * kom2 ) * kom2 ) +
0082                x[ ij ] * ( pow ( ( k / k0[ ij ] ) , l21 ) ) *
0083                G0[ ij ] * ss0[ ij ] /( s0 - ss ) );
0084
0085     fFormula = tangent / sqrt( 1. + pow( tangent , 2 ) );
0086
0087   }
0088
0089   /* !     Note on quadrant corrections: */
0090   /* !     The phaseshifts (delta_Ltj) are from 0-180 degrees. */
0091   /* !     The tangent function is positive from 0-90 degrees and negative from 90-180 deg. */
0092   /* !     Sine is positive from 0-180 deg while Cosine is negative from 90-180 deg. */
0093   /* !     This means that "formula" will give the wrong sign (negative) for sin(delta) */
0094   /* !     from 90-180 deg.  However, this works out okay, since the cosine calculated  */
0095   /* !     in the calling routine is always positive, and the real part of the partial  */
0096   /* !     wave amplitude (cos*sin) will have the correct sign for 0-180 deg. */
0097   /* !     The imaginary part (sin**2) is always positive.   */
0098   /* !     Thus, everything works out correctly. */
0099
0100   return fFormula;
0101
0102 }
0103
0104
0105 /* ! this routine returns the CM angular distribution parameters for pi-N */
0106 /* ! elastic scattering. */
0107
0108 /* ! based on work by A. Shinozaki */
0109 /* ! gh 02.06.20 */
0110
0111 /* ! ich = Isospin channel (1: pure T=3/2,  2: mixed T=1/2,3/2) */
0112 /* ! q  = pi-N CM momentum (GeV/c) */
0113 /* ! s  = pi-N CM energy^2 (GeV) */
0114 /* ! z0 = cos^0 parameter */
0115 /* ! z1 = cos^1 parameter (in pi-N CM frame) */
0116 /* ! z2 = cos^2 parameter */
0117
0118 void phaseshifts( int chan, double q, double s ) {
0119
0120   double hbarc2 = 0.389379; //  mb(GeV/c)^2
0121   double reck2 = hbarc2 / pow( q , 2 );
0122
0123   double ss31 = formula( s , q , 0 , 2 ); // S31
0124   double cs31 = sqrt( 1.0 - pow( ss31 , 2 ) );
0125   std::complex<double> fs31 ( ss31 * cs31 , pow( ss31 , 2 ) );
0126
0127   double sp33 = formula( s , q , 1 , 6 ); // P33
0128   double cp33 = sqrt( 1.0 - pow( sp33 , 2 ) );
0129   std::complex<double> fp33 ( sp33 * cp33 , pow( sp33 , 2 ) );
0130
0131   double sp31 = formula( s , q , 1 , 5 ); // P31
0132   double cp31 = sqrt( 1.0 - pow( sp31 , 2 ) );
0133   std::complex<double> fp31 ( sp31 * cp31 , pow( sp31 , 2 ) );
0134
0135
0136   if ( chan == 1 ) {
0137     std::complex<double> fs11 ( 0., 0. );
0138     std::complex<double> fp11 ( 0., 0. );
0139     std::complex<double> fp13 ( 0., 0. );
0140
0141     // ! equations for A,B,C from Appendix F of Aki's thesis.
0142     // cos^0
0143     fZ0 = ( reck2 *
0144         std::real( ( fs31 + 2.0 * fs11 ) * std::conj( fs31 + 2.0 * fs11 ) +
0145                ( ( fp33 + 2.0 * fp13 ) - ( fp31 + 2.0 * fp11 ) ) *
0146                std::conj( ( fp33 + 2.0 * fp13 ) -
0147                   ( fp31 + 2.0 * fp11 ) ) ) );
0148
0149     // cos^1
0150     fZ1 = ( reck2 *
0151         2.0 * std::real( ( fs31 + 2.0 * fs11 ) *
0152                  std::conj( 2.0 * ( fp33 + 2.0 * fp13 ) +
0153                     ( fp31 + 2.0 * fp11 ) ) ) );
0154
0155
0156     // cos^2
0157     fZ2 = ( reck2 *
0158         3.0 * std::real( ( fp33 + 2.0 * fp13 ) *
0159                  std::conj( fp33 + 2.0 * fp13 ) +
0160                  2.0 * ( fp33 + 2.0 * fp13 ) *
0161                  std::conj( fp31 + 2.0 * fp11 ) ) );
0162
0163   }
0164
0165
0166   // For Pion- . mixed isospin 1/2,3/2 channel (e.g. pi-p)
0167   if ( chan == 2 ) {
0168     double ss11 = formula( s , q , 0 , 1 ); // S11
0169     double cs11 = sqrt( 1.0 - pow( ss11 , 2 ) );
0170     std::complex<double> fs11 ( ss11 * cs11 , pow( ss11 , 2 ) );
0171
0172     double sp13 = formula( s , q , 1 , 4 ); // P13
0173     double cp13 = sqrt( 1.0 - pow( sp13 , 2 ) );
0174     std::complex<double> fp13 ( sp13 * cp13 , pow( sp13  , 2 ) );
0175
0176     double sp11 = formula( s , q , 1 , 3 ); // P11
0177     double cp11 = sqrt( 1.0 - pow( sp11 , 2 ) );
0178     std::complex<double> fp11 ( sp11 * cp11 , pow( sp11 , 2 ) );
0179
0180
0181     // ! equations for A,B,C from Appendix F of Aki's thesis.
0182     // cos^0
0183     fZ0 = ( reck2 *
0184         std::real( ( fs31 + 2.0 * fs11 ) * std::conj( fs31 + 2.0 * fs11 ) +
0185                ( ( fp33 + 2.0 * fp13 ) - ( fp31 + 2.0 * fp11 ) ) *
0186                std::conj( ( fp33 + 2.0 * fp13 ) -
0187                   ( fp31 + 2.0 * fp11 ) ) ) );
0188
0189     // cos^1
0190     fZ1 = ( reck2 *
0191         2.0 * std::real( ( fs31 + 2.0 * fs11 ) *
0192                  std::conj( 2.0 * ( fp33 + 2.0 * fp13 ) +
0193                     ( fp31 + 2.0 * fp11 ) ) ) );
0194
0195
0196     // cos^2
0197     fZ2 = ( reck2 *
0198         3.0 * std::real( ( fp33 + 2.0 * fp13 ) *
0199                  std::conj( fp33 + 2.0 * fp13 ) +
0200                  2.0 * ( fp33 + 2.0 * fp13 ) *
0201                  std::conj( fp31 + 2.0 * fp11 ) ) );
0202
0203     fZ0 = fZ0 / 9;
0204     fZ1 = fZ1 / 9;
0205     fZ2 = fZ2 / 9;
0206   }
0207
0208 }
0209
0210 double getZ0(){ return fZ0; }
0211 double getZ1(){ return fZ1; }
0212 double getZ2(){ return fZ2; }