--- /dev/null
+///////////////////////////////////////////////////////////////////////////
+//
+// Copyright 2010
+//
+// This file is part of starlight.
+//
+// starlight is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// starlight is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with starlight. If not, see <http://www.gnu.org/licenses/>.
+//
+///////////////////////////////////////////////////////////////////////////
+//
+// File and Version Information:
+// $Rev:: 164 $: revision of last commit
+// $Author:: odjuvsla $: author of last commit
+// $Date:: 2013-10-06 16:18:08 +0200 #$: date of last commit
+//
+// Description:
+//
+//
+//
+///////////////////////////////////////////////////////////////////////////
+
+
+#include <iostream>
+#include <fstream>
+
+#include "starlightconstants.h"
+#include "reportingUtils.h"
+#include "nucleus.h"
+#include <inputParameters.h>
+
+
+using namespace std;
+using namespace starlightConstants;
+
+//______________________________________________________________________________
+nucleus::nucleus(const int Z,
+ const int A,
+ const double deuteronSlopePar,
+ const bool dAuCoherentProduction)
+ : _Z(Z),
+ _A(A),
+ _deuteronSlopePar(deuteronSlopePar),
+ _dAuCoherentProduction(dAuCoherentProduction)
+{
+ init();
+}
+
+void nucleus::init()
+{
+ switch (_Z) {
+ case 79:
+ {
+ _Q0 = 0.060;
+ _rho0 = 0.159407;
+ }
+ break;
+ case 53:
+ case 49:
+ {
+ _Q0 = 0.069;
+ _rho0 = 0.161626;
+ }
+ break;
+ case 29:
+ {
+ _Q0 = 0.087;
+ _rho0 = 0.166878;
+ }
+ break;
+ case 14:
+ {
+ _Q0 = 0.115;
+ _rho0 = 0.177128;
+ }
+ break;
+ case 8:
+ {
+ _Q0 = 0.138;
+ _rho0 = 0.188459;
+ }
+ break;
+ case 82:
+ {
+ _Q0 = 0.059;
+ _rho0 = 0.159176;
+ }
+ break;
+ case 20:
+ {
+ _Q0 = 0.102;
+ _rho0 = 0.171907;
+ }
+ break;
+ case 1: // _Q0 and _rho0 are not relevant for protons.
+ {
+ _Q0 = -1.0;
+ _rho0 = -1.0;
+ }
+ break;
+ default:
+ printWarn << "density not defined for projectile with Z = " << _Z << ". using defaults." << endl;
+ _rho0 = 0.16; //'typical' density
+ _Q0 = 0.0; // not used.
+ }
+ _r0 = 1.16 * (1. - 1.16 * pow(_A, -2. / 3.)); // for FRITIOF and FormFactor.
+}
+
+//______________________________________________________________________________
+nucleus::~nucleus()
+{ }
+
+
+//______________________________________________________________________________
+double
+nucleus::nuclearRadius() const
+{
+ // we use this for specific nuclei, Au, Pb, protons...
+ if (_Z == 79) // Au
+ return 6.38; // [fm]
+ if (_Z == 82) // Pb
+ return 6.62; // [fm]
+ if ((_Z == 1) && (_A == 1)) // proton
+ return 0.7; // [fm]
+ if ((_Z == 1) && (_A == 2)) // deuteron
+ return 1.9; // [fm]
+ return 1.2 * pow(_A, 1. / 3.);
+}
+
+
+//______________________________________________________________________________
+double
+nucleus::formFactor(const double t) const
+{
+ // electromagnetic form factor of proton
+ if ((_Z == 1) && (_A == 1)) {
+ const double rec = 1. / (1. + t / 0.71);
+ return rec * rec;
+ }
+ // deuteron form factor
+ if ((_Z == 1) && (_A == 2)) { // careful with this line on dAu
+ // this is for dAu//Sergey
+ // sergey's stuff, also replaced b with _deuteronSlopePar and dropped t02 since it wasnt used
+ // incoherent form factor F(t) = 0.34 e(141.5 t) + 0.58 e(26.1 t) + 0.08 e(15.5 t)
+ const double st = 0.34 * exp(-141.5 * t ) + 0.58 * exp(-26.1 * t ) + 0.08 * exp(-15.5 * t );
+ const double st4 = 0.34 * exp(-141.5 * t / 4) + 0.58 * exp(-26.1 * t / 4) + 0.08 * exp(-15.5 * t / 4);
+ // st paramters from Franco and Varma for st eqn PRL33 ...
+ // form factor from Eisenberg, nuclear physics B 104
+ const double arg = _deuteronSlopePar * t;
+ if (_dAuCoherentProduction)
+ return (st4 * st4 * exp(-arg) - 0.068 * st4 * exp(-arg * 3. / 4.));
+ return exp(-arg) * 0.5 * (1 + st) - 0.068 * exp(-arg * 3. / 4.)
+ - st4 * st4 * exp(-arg) + 0.068 * st4 * exp(-arg * 3. / 4.);
+ }
+ // nuclear form factor
+ // use parameterization from FRITIOF
+ // R = r0 * A^{1/3} with r0 = 1.16 * (1 - 1.16 * A^{-2/3})
+ const double R = fritiofR0();
+ const double q = sqrt(t);
+ const double arg1 = q * R / hbarc;
+ const double arg2 = hbarc / (q * _r0);
+ const double sph = (sin(arg1) - arg1 * cos(arg1)) * 3. * arg2 * arg2 * arg2 / double(_A);
+ const double a0 = 0.70; // [fm]
+ return sph / (1. + (a0 * a0 * t) / (hbarc * hbarc));
+}
+
+//______________________________________________________________________________
+
+double
+nucleus::dipoleFormFactor(const double t, const double t0) const
+{
+ const double rec = 1. / (1. + t / t0);
+ return rec * rec;
+}
+
+//______________________________________________________________________________
+double
+nucleus::thickness(const double b) const
+{
+ // JS This code calculates the nuclear thickness function as per Eq. 4 in
+ // Klein and Nystrand, PRC 60.
+ // former DOUBLE PRECISION FUNCTION T(b)
+
+ // data for Gauss integration
+ const unsigned int nmbPoints = 5;
+ const double xg[nmbPoints + 1] = {0., 0.1488743390, 0.4333953941, 0.6794095683,
+ 0.8650633667, 0.9739065285};
+ const double ag[nmbPoints + 1] = {0., 0.2955242247, 0.2692667193, 0.2190863625,
+ 0.1494513492, 0.0666713443};
+
+ const double zMin = 0;
+ const double zMax = 15;
+ const double zRange = 0.5 * (zMax - zMin);
+ const double zMean = 0.5 * (zMax + zMin);
+ double sum = 0;
+ for(unsigned int i = 1; i <= nmbPoints; ++i) {
+ double zsp = zRange * xg[i] + zMean;
+ double radius = sqrt(b * b + zsp * zsp);
+ sum += ag[i] * rws(radius);
+ zsp = zRange * (-xg[i]) + zMean;
+ radius = sqrt(b * b + zsp * zsp);
+ sum += ag[i] * rws(radius);
+ }
+
+ return 2. * zRange * sum;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff