estarlight/src/nucleus.cpp

0001 ///////////////////////////////////////////////////////////////////////////
0002 //
0003 //    Copyright 2010
0004 //
0005 //    This file is part of starlight.
0006 //
0007 //    starlight is free software: you can redistribute it and/or modify
0008 //    it under the terms of the GNU General Public License as published by
0009 //    the Free Software Foundation, either version 3 of the License, or
0010 //    (at your option) any later version.
0011 //
0012 //    starlight is distributed in the hope that it will be useful,
0013 //    but WITHOUT ANY WARRANTY; without even the implied warranty of
0014 //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0015 //    GNU General Public License for more details.
0016 //
0017 //    You should have received a copy of the GNU General Public License
0018 //    along with starlight. If not, see <http://www.gnu.org/licenses/>.
0019 //
0020 ///////////////////////////////////////////////////////////////////////////
0021 //
0022 // File and Version Information:
0023 // $Rev:: 262                         $: revision of last commit
0024 // $Author:: jnystrand                $: author of last commit
0025 // $Date:: 2016-06-01 14:14:20 +0100 #$: date of last commit
0026 //
0027 // Description:
0028 //
0029 //
0030 //
0031 ///////////////////////////////////////////////////////////////////////////
0032
0033
0034 #include <iostream>
0035 #include <fstream>
0036
0037 #include "starlightconstants.h"
0038 #include "reportingUtils.h"
0039 #include "nucleus.h"
0040 #include <inputParameters.h>
0041
0042
0043 using namespace std;
0044 using namespace starlightConstants;
0045
0046 //______________________________________________________________________________
0047 nucleus::nucleus(const int    Z,
0048                  const int    A,
0049          const int    productionMode)
0050     : _Z(Z),
0051       _A(A),
0052       _productionMode(productionMode)
0053 {
0054   init();
0055 }
0056
0057 void nucleus::init()
0058 {
0059   switch (_Z) {
0060     case 82:
0061         {
0062           _Radius = 6.624;
0063           _rho0 = 0.160696;
0064         }
0065         break;
0066     case 79:
0067         {
0068           _Radius = 6.38;
0069           _rho0 = 0.169551;
0070         }
0071         break;
0072     case 29:
0073         {
0074                   _Radius = 4.214;
0075           _rho0 = 0.173845;
0076         }
0077         break;
0078     case 1:
0079         {
0080           //is this a proton or deuteron
0081           if(_A==1){
0082             _Radius = 0.7;
0083             _rho0 = -1.0; //Not relevant for protons
0084           }
0085           // this is electron
0086           else if(_A==0){
0087             _Radius = 0.;
0088             _rho0 = -1.0;
0089           }
0090           else {
0091             _Radius = 2.1;
0092             _rho0 = _A;
0093           }
0094         }
0095         break;
0096     case -1:
0097         {
0098           //is this a anti-proton
0099           if(_A==1){
0100             _Radius = 0.7;
0101             _rho0 = -1.0; //Not relevant for protons
0102           }
0103           // this is positron
0104           else if(_A==0){
0105             _Radius = 0.;
0106             _rho0 = -1.0;
0107           }
0108           else {
0109             _Radius = 2.1;
0110             _rho0 = _A;
0111           }
0112         }
0113         break;
0114     default:
0115         printWarn << "density not defined for projectile with Z = " << _Z << ". using defaults." << endl;
0116                 _Radius = 1.2*pow(_A, 1. / 3.);
0117         _rho0 = 0.138/(1.13505-0.0004283*_A);  //This matches the radius above
0118         if( _Z < 7 ){
0119           // This is for Gaussian form factors/densities
0120           _rho0 = _A;
0121         }
0122     }
0123 }
0124
0125 //______________________________________________________________________________
0126 nucleus::~nucleus()
0127 { }
0128
0129 //______________________________________________________________________________
0130 double
0131 nucleus::rws(const double r) const
0132 {
0133   if( _Z < 7 ){
0134     // Gaussian density distribution for light nuclei
0135     double norm = (3./(2.*starlightConstants::pi))*sqrt( (3./(2.*starlightConstants::pi)) );
0136     norm = norm/(nuclearRadius()*nuclearRadius()*nuclearRadius());
0137     return norm*exp(-((3./2.)*r*r)/(nuclearRadius()*nuclearRadius()));
0138   }else{
0139     // Fermi density distribution for heavy nuclei
0140     return 1.0 / (1. + exp((r - nuclearRadius()) / woodSaxonSkinDepth()));
0141   }
0142 }
0143
0144 //______________________________________________________________________________
0145 double
0146 nucleus::formFactor(const double t) const
0147 {
0148     // electromagnetic form factor of proton
0149     if ((_Z == 1) && (_A == 1)) {
0150         const double rec = 1. / (1. + t / 0.71);
0151         return rec * rec;
0152     }
0153
0154         if( _Z < 7 ){
0155       // Gaussian form factor for light nuclei
0156           const double R_G    = nuclearRadius();
0157           return exp(-R_G*R_G*t/(6.*starlightConstants::hbarc*starlightConstants::hbarc));
0158         } else {
0159           // nuclear form factor, from Klein Nystrand PRC 60 (1999) 014903, Eq. 14
0160           const double R    = nuclearRadius();
0161           const double q    = sqrt(t);
0162           const double arg1 = q * R / hbarc;
0163           const double arg2 = hbarc / (q * R);
0164           const double sph  = (sin(arg1) - arg1 * cos(arg1)) * 3. * arg2 * arg2 * arg2;
0165           const double a0   = 0.70;  // [fm]
0166           return sph / (1. + (a0 * a0 * t) / (hbarc * hbarc));
0167         }
0168 }
0169
0170 //______________________________________________________________________________
0171
0172 double
0173 nucleus::dipoleFormFactor(const double t, const double t0) const
0174 {
0175      const double rec = 1. / (1. + t / t0);
0176      return rec * rec;
0177 }
0178
0179 //______________________________________________________________________________
0180 double
0181 nucleus::thickness(const double b) const
0182 {
0183     //    JS      This code calculates the nuclear thickness function as per Eq. 4 in
0184     //    Klein and Nystrand, PRC 60.
0185     //    former DOUBLE PRECISION FUNCTION T(b)
0186
0187     // data for Gauss integration
0188     const unsigned int nmbPoints         = 5;
0189     const double       xg[nmbPoints + 1] = {0., 0.1488743390, 0.4333953941, 0.6794095683,
0190                                             0.8650633667, 0.9739065285};
0191     const double       ag[nmbPoints + 1] = {0., 0.2955242247, 0.2692667193, 0.2190863625,
0192                                             0.1494513492, 0.0666713443};
0193
0194     const double zMin   = 0;
0195     const double zMax   = 15;
0196     const double zRange = 0.5 * (zMax - zMin);
0197     const double zMean  = 0.5 * (zMax + zMin);
0198     double       sum    = 0;
0199     for(unsigned int i = 1; i <= nmbPoints; ++i) {
0200         double zsp    = zRange * xg[i] + zMean;
0201         double radius = sqrt(b * b + zsp * zsp);
0202         sum          += ag[i] * rws(radius);
0203         zsp           = zRange * (-xg[i]) + zMean;
0204         radius        = sqrt(b * b + zsp * zsp);
0205         sum          += ag[i] * rws(radius);
0206     }
0207
0208     return 2. * zRange * sum;
0209 }