--- /dev/null
+
+/*
+ <one line to give the program's name and a brief idea of what it does.>
+ Copyright (C) <year> <name of author>
+
+p This program 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.
+
+ This program 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 this program. If not, see <http://www.gnu.org/licenses/>.
+
+*/
+
+#include "spectrum.h"
+#include <cmath>
+#include "beambeamsystem.h"
+#include <randomgenerator.h>
+#include <iostream>
+
+spectrum::spectrum(beamBeamSystem *bbs) :
+ _bMin(5.0)
+ ,_bMax(128000.0)
+ ,_nBbins(6400)
+ ,_probOfBreakup(_nBbins)
+ ,_beamBeamSystem(bbs)
+ ,_nK(10000)
+ ,_fnSingle(_nK)
+ ,_fnDouble(_nK)
+ ,_fnSingleCumulative(_nK+1)
+ ,_fnDoubleCumulative(_nK+1)
+ ,_fnDoubleInt(_nK)
+ ,_fnDoubleIntCumulative(_nK+1)
+ ,_eGamma(_nK+1)
+ ,_eGammaMin(6.0)
+ ,_eGammaMax(600000.0)
+ ,_zTarget(82)
+ ,_aTarget(278)
+ ,_hadBreakProbCalculated(false)
+{
+ _eGamma.resize(_nK+1);
+ _probOfBreakup.resize(_nBbins);
+}
+
+int spectrum::generateKsingle()
+{
+
+ _fnSingle.resize(_nK);
+ _fnSingleCumulative.resize(_nK+1);
+
+ double eg_inc = exp(log(_eGammaMax/_eGammaMin)/(double)_nK);
+
+ double egamma = _eGammaMin;
+ for (int i = 0; i < _nK+1; i++)
+ {
+ _eGamma[i] = egamma;
+ egamma = egamma * eg_inc;
+ }
+ egamma = _eGammaMin;
+
+ double fnorm = 0;
+
+
+ if (_hadBreakProbCalculated == false)
+ {
+ _hadBreakProbCalculated = generateBreakupProbabilities();
+ }
+ double binc = exp((log(_bMax/_bMin))/(double)_nBbins);
+
+ for (int i = 0; i < _nK; i++)
+ {
+ double b = _bMin;
+
+ double bint = 0.0;
+
+ double f1 = 0;
+ double f2 = 0;
+
+ for (int j = 0; j < _nBbins - 1; j++)
+ {
+ double bold = b;
+ if (j == 0)
+ {
+ //f1 = fBeamBeamSystem->getBeam1().nofe(egamma, b)*GetSigma(egamma)*fProbOfBreakup[j]*b;
+ f1 = getTransformedNofe(egamma, b)*getSigma(egamma)*_probOfBreakup[j]*b;
+ //std::cout << fProbOfBreakup[j] << std::endl;
+ }
+ else
+ {
+ f1 = f2;
+ }
+ b = b*binc;
+// f2 = fBeamBeamSystem->getBeam1().nofe(egamma, b)*GetSigma(egamma)*fProbOfBreakup[j+1]*b;;
+ f2 = getTransformedNofe(egamma, b)*getSigma(egamma)*_probOfBreakup[j+1]*b;;
+ bint = bint + 0.5*(f1+f2)*(b-bold);
+ }
+ bint = 2.0*starlightConstants::pi*bint;
+ if (i == 0)
+ {
+ fnorm = 1.0/bint;
+ }
+ _fnSingle[i] = bint*(_eGamma[i+1]-_eGamma[i]);
+
+ egamma = egamma*eg_inc;
+ }
+
+ _fnSingleCumulative[0] = 0.00;
+ for (int i = 0; i < _nK; i++)
+ {
+ _fnSingleCumulative[i+1] = _fnSingleCumulative[i]+_fnSingle[i];
+ }
+
+ double fnormfactor = 1.00/_fnSingleCumulative[_nK];
+ for (int i = 0; i < _nK; i++)
+ {
+ _fnSingleCumulative[i+1] = fnormfactor*_fnSingleCumulative[i+1];
+ }
+
+ return 0;
+
+}
+
+int spectrum::generateKdouble()
+{
+ //Quick fix for now TODO: Fix it!
+ _nK = 100;
+
+ _fnDouble.resize(_nK);
+ _fnDoubleInt.resize(_nK);
+ _fnDoubleIntCumulative.resize(_nK+1);
+ _fnDoubleCumulative.resize(_nK+1);
+ for (int i = 0; i < _nK; i++)
+ {
+ _fnDouble[i].resize(_nK);
+ _fnDoubleCumulative[i].resize(_nK+1);
+ }
+ _fnDoubleCumulative[_nK].resize(_nK+1);
+
+ double eg_inc = exp(log(_eGammaMax/_eGammaMin)/(double)_nK);
+ double egamma1 = _eGammaMin;
+ double egamma2 = _eGammaMin;
+
+ for (int i = 0; i < _nK+1; i++)
+ {
+ _eGamma[i] = egamma1;
+ egamma1 = egamma1 * eg_inc;
+ }
+ egamma1 = _eGammaMin;
+
+ double fnorm = 0;
+
+ if (_hadBreakProbCalculated == false)
+ {
+ _hadBreakProbCalculated = generateBreakupProbabilities();
+ }
+
+ double binc = exp((log(_bMax/_bMin))/(double)_nBbins);
+
+ int nbbins = _nBbins;
+
+ //double b_min = _bMin;
+ //double b_max = _bMax;
+
+ for (int i = 0; i < _nK; i++)
+ {
+
+ egamma2 = _eGammaMin;
+ //double sum_over_k = 0.0;
+
+ for (int j = 0; j < _nK; j++)
+ {
+ double bint = 0.0;
+ double b = _bMin;
+ double f1 = 0;
+ double f2 = 0;
+
+ for (int k = 0; k < nbbins - 1; k++)
+ {
+ double bold = b;
+
+ if (k == 0)
+ {
+ // f1 = fBeamBeamSystem->getBeam1().nofe(egamma1, b) * fBeamBeamSystem->getBeam2().nofe(egamma2, b)
+ // * GetSigma(egamma1) * GetSigma(egamma2) *fProbOfBreakup[j]*b;
+ f1 = getTransformedNofe(egamma1, b) * getTransformedNofe(egamma2, b)
+ * getSigma(egamma1) * getSigma(egamma2) *_probOfBreakup[k]*b; }
+ else
+ {
+ f1 = f2;
+ }
+ b = b*binc;
+ // f2 = fBeamBeamSystem->getBeam1().nofe(egamma1, b) * fBeamBeamSystem->getBeam2().nofe(egamma2, b)
+ // * GetSigma(egamma1) * GetSigma(egamma2) *fProbOfBreakup[j+1]*b;
+ f2 = getTransformedNofe(egamma1, b) * getTransformedNofe(egamma2, b)
+ * getSigma(egamma1) * getSigma(egamma2) *_probOfBreakup[k+1]*b;
+ bint = bint + 0.5*(f1+f2)*(b-bold);
+ }
+ bint = 2.0*starlightConstants::pi*bint;
+ _fnDouble[i][j] = bint * (_eGamma[i+1] - _eGamma[i]) * (_eGamma[j+1] - _eGamma[j]);
+ egamma2 = egamma2 * eg_inc;
+ }
+ egamma1 = egamma1 * eg_inc;
+ }
+
+ for (int i = 0; i < _nK; i++)
+ {
+ _fnDoubleInt[i] = 0.0;
+ for (int j = 0; j < _nK; j++)
+ {
+ _fnDoubleInt[i] = _fnDoubleInt[i] + _fnDouble[i][j];
+ }
+ }
+
+ _fnDoubleIntCumulative[0] = 0.0;
+ for (int i = 1; i < _nK+1; i++)
+ {
+ _fnDoubleIntCumulative[i] = _fnDoubleIntCumulative[i-1] + _fnDoubleInt[i-1];
+ }
+
+ fnorm = 1.0/_fnDoubleIntCumulative[_nK];
+ for (int i = 0; i < _nK+1; i++)
+ {
+ _fnDoubleIntCumulative[i] = fnorm * _fnDoubleIntCumulative[i];
+ }
+
+ return 0;
+}
+
+double spectrum::drawKsingle()
+{
+ double xtest = 0;
+ int itest = 0;
+ double egamma = 0.0;
+
+ xtest = randyInstance.Rndom();
+ while (xtest > _fnSingleCumulative[itest])
+ {
+ itest++;
+ }
+ itest = itest - 1;
+
+ if (itest >= _nK || itest < 0)
+ {
+ std::cerr << "ERROR: itest: " << itest << std::endl;
+
+ }
+
+ double delta_f = xtest - _fnSingleCumulative[itest];
+ if (delta_f <= 0.0)
+ {
+ std::cout << "WARNING: delta_f: " << delta_f << std::endl;
+ return -1;
+ }
+ double dE = _eGamma[itest+1] - _eGamma[itest];
+ double dF = _fnSingleCumulative[itest+1] - _fnSingleCumulative[itest];
+
+ double delta_e = delta_f*dE/dF;
+
+ if (delta_e > (_eGamma[itest+1] - _eGamma[itest]))
+ {
+ std::cerr << "ERROR: delta_E: " << delta_e << std::endl;
+ }
+
+ egamma = _eGamma[itest] + delta_e;
+ return egamma;
+}
+
+void spectrum::drawKdouble(float& egamma1, float& egamma2)
+{
+ double xtest1 = 0.0;
+ double xtest2 = 0.0;
+ int itest1 = 0;
+ int itest2 = 0;
+
+ xtest1 = randyInstance.Rndom();
+
+ while (xtest1 > _fnDoubleIntCumulative[itest1])
+ {
+ itest1++;
+ }
+ itest1 = itest1 - 1;
+
+ if (itest1 >= _nK || itest1 < 0)
+ {
+ std::cerr << "ERROR: itest1: " << itest1 << std::endl;
+ }
+ double delta_f = xtest1 - _fnDoubleIntCumulative[itest1];
+
+ if (delta_f <= 0.0)
+ {
+ std::cout << "WARNING: delta_f: " << delta_f << std::endl;
+ }
+
+ double dE = _eGamma[itest1+1] - _eGamma[itest1];
+ double dF = _fnDoubleIntCumulative[itest1+1] - _fnDoubleIntCumulative[itest1];
+
+ double delta_e = delta_f*dE/dF;
+
+ if (delta_e > (_eGamma[itest1+1] - _eGamma[itest1]))
+ {
+ std::cerr << "ERROR: delta_E: " << delta_e << std::endl;
+ }
+
+ egamma1 = _eGamma[itest1] + delta_e;
+
+ // Second gamma
+
+ // double reldw = delta_e/(fEGamma[itest1+1] - fEGamma[itest1]);
+ // std::vector<double> fn_second(fNK);
+ std::vector<double> fn_second_cumulative(_nK+1);
+
+ // for(int i = 0; i < fNK; i++)
+ // {
+ // fn_second[i] = fFnDouble[itest1][i] + (fFnDouble[itest1+1][i] - fFnDouble[itest1][i])*reldw;
+ // }
+
+ fn_second_cumulative[0] = 0.0;
+ for(int i = 1; i < _nK+1; i++)
+ {
+ // fn_second_cumulative[i] = fn_second_cumulative[i-1] + fn_second[i-1]; //TODO:check indexing
+ fn_second_cumulative[i] = fn_second_cumulative[i-1] + _fnDouble[itest1][i-1];
+ }
+
+ double norm_factor = 1.0/fn_second_cumulative[_nK];
+ for(int i = 0; i < _nK+1; i++)
+ {
+ fn_second_cumulative[i] = norm_factor*fn_second_cumulative[i];
+ }
+
+ xtest2 = randyInstance.Rndom();
+
+ while (xtest2 > fn_second_cumulative[itest2])
+ {
+ itest2++;
+ }
+ itest2 = itest2 - 1;
+
+ if (itest2 >= _nK || itest2 < 0)
+ {
+ std::cerr << "ERROR: itest2: " << itest2 << std::endl;
+ }
+ delta_f = xtest2 - fn_second_cumulative[itest2];
+
+ if (delta_f <= 0.0)
+ {
+ std::cout << "WARNING: delta_f: " << delta_f << std::endl;
+ }
+
+ dE = _eGamma[itest2+1] - _eGamma[itest2];
+ dF = fn_second_cumulative[itest2+1] - fn_second_cumulative[itest2];
+
+ delta_e = delta_f*dE/dF;
+
+ if (delta_e > (_eGamma[itest2+1] - _eGamma[itest2]))
+ {
+ std::cerr << "ERROR: delta_E: " << delta_e << std::endl;
+ }
+
+ egamma2 = _eGamma[itest2] + delta_e;
+
+ return;
+}
+
+
+bool spectrum::generateBreakupProbabilities()
+{
+
+ int nbbins = _nBbins;
+
+ double b_min = _bMin;
+ //double b_max = _bMax;
+
+ // double binc = (log(b_max/b_min))/(double)nbbins;
+ double binc = exp((log(_bMax/_bMin))/(double)_nBbins);
+
+ double b = b_min;
+
+ _probOfBreakup.resize(nbbins);
+
+ for (int i = 0; i < nbbins; i++)
+ {
+ double bimp = b;
+ double rhad = 0;
+ rhad = _beamBeamSystem->probabilityOfBreakup(bimp);
+ _probOfBreakup[i] = exp(-rhad);
+ b = b*binc;
+ }
+ return true;
+}
+
+double spectrum::getFnSingle(double egamma) const
+{
+ double eginc = exp(log(_eGammaMax/_eGammaMin)/static_cast<double>(_nK));
+ int i1 = log(egamma/_eGammaMin)/log(eginc);
+ int i2 = i1 + 1;
+ double fnSingle = 0.0;
+
+ if (i1 < 0 || i2 > _nK)
+ {
+ std::cout << "I1, I2 out of bounds. Egamma = " << egamma << std::endl;
+ std::cout << "I1, I2 = " << i1 << ", " << i2 << std::endl;
+ fnSingle = 0.0;
+ }
+ else
+ {
+ double dE = _eGamma[i2] - _eGamma[i1];
+ double eFrac = (egamma - _eGamma[i1])/dE;
+
+ if (eFrac < 0.0 || eFrac > 1.0)
+ {
+ std::cout << "WARNING: Efrac = " << eFrac << std::endl;
+ }
+ fnSingle = (1.0 - eFrac)*_fnSingle[i1] + eFrac*_fnSingle[i2];
+ }
+ return fnSingle;
+}
+
+double spectrum::getFnDouble(double egamma1, double egamma2) const
+{
+ double eginc = exp(log(_eGammaMax/_eGammaMin)/static_cast<double>(_nK));
+ int i1 = log(egamma1/_eGammaMin)/log(eginc);
+ int i2 = i1 + 1;
+ double fnDouble = 0.0;
+
+ if (i1 < 0 || i2 > _nK)
+ {
+ std::cout << "I1, I2 out of bounds. Egamma1 = " << egamma1 << std::endl;
+ std::cout << "I1, I2 = " << i1 << ", " << i2 << std::endl;
+ fnDouble = 0.0;
+ return fnDouble;
+ }
+
+ int j1 = log(egamma2/_eGammaMin)/log(eginc);
+ int j2 = j1 + 1;
+
+ if (j1 < 0 || j2 > _nK)
+ {
+ std::cout << "J1, J2 out of bounds. Egamma2 = " << egamma2 << std::endl;
+ std::cout << "J1, J2 = " << j1 << ", " << j2 << std::endl;
+ fnDouble = 0.0;
+ return fnDouble;
+ }
+
+ double dE1 = _eGamma[i2] - _eGamma[i1];
+ double eFrac1 = (egamma1 - _eGamma[i1])/dE1;
+
+ if (eFrac1 < 0.0 || eFrac1 > 1.0)
+ {
+ std::cout << "WARNING: Efrac1 = " << eFrac1 << std::endl;
+ }
+
+ double ptemp1 = (1.0 - eFrac1)*_fnDouble[i1][j1] + eFrac1*_fnDouble[i2][j1];
+ double ptemp2 = (1.0 - eFrac1)*_fnDouble[i1][j2] + eFrac1*_fnDouble[i2][j2];
+
+ double dE2 = _eGamma[j2] - _eGamma[j1];
+ double eFrac2 = (egamma2 - _eGamma[j1])/dE2;
+
+ if (eFrac2 < 0.0 || eFrac2 > 1.0)
+ {
+ std::cout << "WARNING: Efrac2 = " << eFrac2 << std::endl;
+ }
+
+ fnDouble = (1.0 - eFrac2)*ptemp1 + eFrac2*ptemp2;
+
+ return fnDouble;
+
+}
+
+double spectrum::getTransformedNofe(double egamma, double b)
+{
+ double factor = 1.0/(2.0*_beamBeamSystem->beamLorentzGamma());
+ double res = factor * _beamBeamSystem->beam1().photonFlux(b, egamma*factor);
+
+ return res;
+}
+
+
+
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff