STARLIGHT/starlight/src/.svn/text-base/nucleus.cpp.svn-base

   1 ///////////////////////////////////////////////////////////////////////////
   2 //
   3 //    Copyright 2010
   4 //
   5 //    This file is part of starlight.
   6 //
   7 //    starlight is free software: you can redistribute it and/or modify
   8 //    it under the terms of the GNU General Public License as published by
   9 //    the Free Software Foundation, either version 3 of the License, or
  10 //    (at your option) any later version.
  11 //
  12 //    starlight is distributed in the hope that it will be useful,
  13 //    but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  15 //    GNU General Public License for more details.
  16 //
  17 //    You should have received a copy of the GNU General Public License
  18 //    along with starlight. If not, see <http://www.gnu.org/licenses/>.
  19 //
  20 ///////////////////////////////////////////////////////////////////////////
  21 //
  22 // File and Version Information:
  23 // $Rev::                             $: revision of last commit
  24 // $Author::                          $: author of last commit
  25 // $Date::                            $: date of last commit
  26 //
  27 // Description:
  28 //
  29 //
  30 //
  31 ///////////////////////////////////////////////////////////////////////////
  32
  33
  34 #include <iostream>
  35 #include <fstream>
  36
  37 #include "starlightconstants.h"
  38 #include "reportingUtils.h"
  39 #include "nucleus.h"
  40 #include <inputParameters.h>
  41
  42
  43 using namespace std;
  44 using namespace starlightConstants;
  45
  46 //______________________________________________________________________________
  47 nucleus::nucleus(const int    Z,
  48                  const int    A,
  49                  const double deuteronSlopePar,
  50                  const bool   dAuCoherentProduction)
  51         : _Z(Z),
  52           _A(A),
  53           _deuteronSlopePar(deuteronSlopePar),
  54           _dAuCoherentProduction(dAuCoherentProduction)
  55 {
  56   init();
  57 }
  58
  59 void nucleus::init()
  60 {
  61   switch (_Z) {
  62         case 79:
  63                 {
  64                         _Q0   = 0.060;
  65                         _rho0 = 0.159407;
  66                 }
  67                 break;
  68         case 53:
  69         case 49:
  70                 {
  71                         _Q0   = 0.069;
  72                         _rho0 = 0.161626;
  73                 }
  74                 break;
  75         case 29:
  76                 {
  77                         _Q0   = 0.087;
  78                         _rho0 = 0.166878;
  79                 }
  80                 break;
  81         case 14:
  82                 {
  83                         _Q0   = 0.115;
  84                         _rho0 = 0.177128;
  85                 }
  86                 break;
  87         case 8:
  88                 {
  89                         _Q0   = 0.138;
  90                         _rho0 = 0.188459;
  91                 }
  92                 break;
  93         case 82:
  94                 {
  95                         _Q0   = 0.059;
  96                         _rho0 = 0.159176;
  97                 }
  98                 break;
  99         case 20:
 100                 {
 101                         _Q0   = 0.102;
 102                         _rho0 = 0.171907;
 103                 }
 104                 break;
 105         case 1: // _Q0 and _rho0 are not relevant for protons.
 106                 {
 107                         _Q0   = -1.0;
 108                         _rho0 = -1.0;
 109                 }
 110                 break;
 111         default:
 112                 printWarn << "density not defined for projectile with Z = " << _Z << ". using defaults." << endl;
 113                 _rho0 = 0.16;  //'typical' density
 114                 _Q0   = 0.0;   // not used.
 115         }
 116         _r0 = 1.16 * (1. - 1.16 * pow(_A, -2. / 3.));  // for FRITIOF and FormFactor.
 117 }
 118
 119 //______________________________________________________________________________
 120 nucleus::~nucleus()
 121 { }
 122
 123
 124 //______________________________________________________________________________
 125 double
 126 nucleus::nuclearRadius() const
 127 {
 128         // we use this for specific nuclei, Au, Pb, protons...
 129         if (_Z == 79)                // Au
 130                 return 6.38;  // [fm]
 131         if (_Z == 82)                // Pb
 132                 return 6.62;  // [fm]
 133         if ((_Z == 1) && (_A == 1))  // proton
 134                 return 0.7;   // [fm]
 135         if ((_Z == 1) && (_A == 2))  // deuteron
 136                 return 1.9;   // [fm]
 137         return 1.2 * pow(_A, 1. / 3.);
 138 }
 139
 140
 141 //______________________________________________________________________________
 142 double
 143 nucleus::formFactor(const double t) const
 144 {
 145         // electromagnetic form factor of proton
 146         if ((_Z == 1) && (_A == 1)) {
 147                 const double rec = 1. / (1. + t / 0.71);
 148                 return rec * rec;
 149         }
 150         // deuteron form factor
 151         if ((_Z == 1) && (_A == 2)) {   // careful with this line on dAu
 152                 // this is for dAu//Sergey
 153                 // sergey's stuff, also replaced b with _deuteronSlopePar and dropped t02 since it wasnt used
 154                 // incoherent form factor F(t) = 0.34 e(141.5 t) + 0.58 e(26.1 t) + 0.08 e(15.5 t)
 155                 const double st  = 0.34 * exp(-141.5 * t    ) + 0.58 * exp(-26.1 * t    ) + 0.08 * exp(-15.5 * t    );
 156                 const double st4 = 0.34 * exp(-141.5 * t / 4) + 0.58 * exp(-26.1 * t / 4) + 0.08 * exp(-15.5 * t / 4);
 157                 // st paramters from Franco and Varma for st eqn PRL33 ...
 158                 // form factor from Eisenberg, nuclear physics B 104
 159                 const double arg = _deuteronSlopePar * t;
 160                 if (_dAuCoherentProduction)
 161                         return (st4 * st4 * exp(-arg) - 0.068 * st4 * exp(-arg * 3. / 4.));
 162                 return exp(-arg) * 0.5 * (1 + st) - 0.068 * exp(-arg * 3. / 4.)
 163                         - st4 * st4 * exp(-arg) + 0.068 * st4 * exp(-arg * 3. / 4.);
 164         }
 165         // nuclear form factor
 166         // use parameterization from FRITIOF
 167         // R = r0 * A^{1/3} with r0 = 1.16 * (1 - 1.16 * A^{-2/3})
 168         const double R    = fritiofR0();
 169         const double q    = sqrt(t);
 170         const double arg1 = q * R / hbarc;
 171         const double arg2 = hbarc / (q * _r0);
 172         const double sph  = (sin(arg1) - arg1 * cos(arg1)) * 3. * arg2 * arg2 * arg2 / double(_A);
 173         const double a0   = 0.70;  // [fm]
 174         return sph / (1. + (a0 * a0 * t) / (hbarc * hbarc));
 175 }
 176
 177 //______________________________________________________________________________
 178
 179 double
 180 nucleus::dipoleFormFactor(const double t, const double t0) const
 181 {
 182      const double rec = 1. / (1. + t / t0);
 183      return rec * rec;
 184 }
 185
 186 //______________________________________________________________________________
 187 double
 188 nucleus::thickness(const double b) const
 189 {
 190         //    JS      This code calculates the nuclear thickness function as per Eq. 4 in
 191         //    Klein and Nystrand, PRC 60.
 192         //    former DOUBLE PRECISION FUNCTION T(b)
 193
 194         // data for Gauss integration
 195         const unsigned int nmbPoints         = 5;
 196         const double       xg[nmbPoints + 1] = {0., 0.1488743390, 0.4333953941, 0.6794095683,
 197                                                 0.8650633667, 0.9739065285};
 198         const double       ag[nmbPoints + 1] = {0., 0.2955242247, 0.2692667193, 0.2190863625,
 199                                                 0.1494513492, 0.0666713443};
 200
 201         const double zMin   = 0;
 202         const double zMax   = 15;
 203         const double zRange = 0.5 * (zMax - zMin);
 204         const double zMean  = 0.5 * (zMax + zMin);
 205         double       sum    = 0;
 206         for(unsigned int i = 1; i <= nmbPoints; ++i) {
 207                 double zsp    = zRange * xg[i] + zMean;
 208                 double radius = sqrt(b * b + zsp * zsp);
 209                 sum          += ag[i] * rws(radius);
 210                 zsp           = zRange * (-xg[i]) + zMean;
 211                 radius        = sqrt(b * b + zsp * zsp);
 212                 sum          += ag[i] * rws(radius);
 213         }
 214
 215         return 2. * zRange * sum;
 216 }