estarlight/src/spectrum.cpp

0001
0002 /*
0003     <one line to give the program's name and a brief idea of what it does.>
0004     Copyright (C) <year>  <name of author>
0005
0006 p    This program is free software: you can redistribute it and/or modify
0007     it under the terms of the GNU General Public License as published by
0008     the Free Software Foundation, either version 3 of the License, or
0009     (at your option) any later version.
0010
0011     This program is distributed in the hope that it will be useful,
0012     but WITHOUT ANY WARRANTY; without even the implied warranty of
0013     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0014     GNU General Public License for more details.
0015
0016     You should have received a copy of the GNU General Public License
0017     along with this program.  If not, see <http://www.gnu.org/licenses/>.
0018
0019 */
0020
0021 #include "spectrum.h"
0022 #include <cmath>
0023 #include "beambeamsystem.h"
0024 #include <randomgenerator.h>
0025 #include <iostream>
0026
0027 spectrum::spectrum(const randomGenerator &randy, beamBeamSystem *bbs) :
0028      _bMin(5.0)
0029         ,_bMax(128000.0)
0030     ,_nBbins(6400)
0031     ,_probOfBreakup(_nBbins)
0032         ,_beamBeamSystem(bbs)
0033     ,_nK(10000)
0034         ,_fnSingle(_nK)
0035         ,_fnDouble(_nK)
0036         ,_fnSingleCumulative(_nK+1)
0037         ,_fnDoubleCumulative(_nK+1)
0038         ,_fnDoubleInt(_nK)
0039         ,_fnDoubleIntCumulative(_nK+1)
0040         ,_eGamma(_nK+1)
0041         ,_eGammaMin(6.0)
0042         ,_eGammaMax(600000.0)
0043         ,_zTarget(82)
0044         ,_aTarget(278)
0045         ,_hadBreakProbCalculated(false)
0046     ,_randy(randy)
0047 {
0048     _eGamma.resize(_nK+1);
0049     _probOfBreakup.resize(_nBbins);
0050 }
0051
0052 int spectrum::generateKsingle()
0053 {
0054
0055     _fnSingle.resize(_nK);
0056     _fnSingleCumulative.resize(_nK+1);
0057
0058     double eg_inc = exp(log(_eGammaMax/_eGammaMin)/(double)_nK);
0059
0060     double egamma =  _eGammaMin;
0061     for (int i = 0; i < _nK+1; i++)
0062     {
0063         _eGamma[i] = egamma;
0064         egamma = egamma * eg_inc;
0065     }
0066     egamma = _eGammaMin;
0067
0068     double fnorm = 0;
0069
0070
0071     if (_hadBreakProbCalculated == false)
0072     {
0073         _hadBreakProbCalculated = generateBreakupProbabilities();
0074     }
0075     double binc = exp((log(_bMax/_bMin))/(double)_nBbins);
0076
0077     for (int i = 0; i < _nK; i++)
0078     {
0079         double b = _bMin;
0080
0081         double bint = 0.0;
0082
0083         double f1 = 0;
0084         double f2 = 0;
0085
0086         for (int j = 0; j < _nBbins - 1; j++)
0087         {
0088             double bold = b;
0089             if (j == 0)
0090             {
0091         f1 = getTransformedNofe(egamma, b)*getSigma(egamma)*_probOfBreakup[j]*b;
0092         //std::cout << fProbOfBreakup[j] << std::endl;
0093             }
0094             else
0095             {
0096                 f1 = f2;
0097             }
0098             b = b*binc;
0099             f2 = getTransformedNofe(egamma, b)*getSigma(egamma)*_probOfBreakup[j+1]*b;;
0100             bint = bint + 0.5*(f1+f2)*(b-bold);
0101         }
0102         bint = 2.0*starlightConstants::pi*bint;
0103         if (i == 0)
0104         {
0105             fnorm = 1.0/bint;
0106         }
0107         _fnSingle[i] = bint*(_eGamma[i+1]-_eGamma[i]);
0108
0109         egamma = egamma*eg_inc;
0110     }
0111
0112     _fnSingleCumulative[0] = 0.00;
0113     for (int i = 0; i < _nK; i++)
0114     {
0115         _fnSingleCumulative[i+1] = _fnSingleCumulative[i]+_fnSingle[i];
0116     }
0117
0118     double fnormfactor = 1.00/_fnSingleCumulative[_nK];
0119     for (int i = 0; i < _nK; i++)
0120     {
0121         _fnSingleCumulative[i+1] = fnormfactor*_fnSingleCumulative[i+1];
0122     }
0123
0124     return 0;
0125
0126 }
0127
0128 int spectrum::generateKdouble()
0129 {
0130     //Quick fix for now TODO: Fix it!
0131     _nK = 100;
0132
0133     _fnDouble.resize(_nK);
0134     _fnDoubleInt.resize(_nK);
0135     _fnDoubleIntCumulative.resize(_nK+1);
0136     _fnDoubleCumulative.resize(_nK+1);
0137     for (int i = 0; i < _nK; i++)
0138     {
0139         _fnDouble[i].resize(_nK);
0140         _fnDoubleCumulative[i].resize(_nK+1);
0141     }
0142     _fnDoubleCumulative[_nK].resize(_nK+1);
0143
0144     double eg_inc = exp(log(_eGammaMax/_eGammaMin)/(double)_nK);
0145     double egamma1 =  _eGammaMin;
0146     double egamma2 =  _eGammaMin;
0147
0148     for (int i = 0; i < _nK+1; i++)
0149     {
0150         _eGamma[i] = egamma1;
0151         egamma1 = egamma1 * eg_inc;
0152     }
0153     egamma1 = _eGammaMin;
0154
0155     double fnorm = 0;
0156
0157     if (_hadBreakProbCalculated == false)
0158     {
0159         _hadBreakProbCalculated = generateBreakupProbabilities();
0160     }
0161
0162     double binc = exp((log(_bMax/_bMin))/(double)_nBbins);
0163
0164     int nbbins = _nBbins;
0165
0166     for (int i = 0; i < _nK; i++)
0167     {
0168
0169         egamma2 = _eGammaMin;
0170
0171         for (int j = 0; j < _nK; j++)
0172         {
0173             double bint = 0.0;
0174             double b = _bMin;
0175             double f1 = 0;
0176             double f2 = 0;
0177
0178             for (int k = 0; k < nbbins - 1; k++)
0179             {
0180                 double bold = b;
0181
0182                 if (k == 0)
0183                 {
0184                   f1 = getTransformedNofe(egamma1, b) * getTransformedNofe(egamma2, b)
0185                          * getSigma(egamma1) * getSigma(egamma2) *_probOfBreakup[k]*b;                }
0186                 else
0187                 {
0188                     f1 = f2;
0189                 }
0190                 b = b*binc;
0191                 f2 = getTransformedNofe(egamma1, b) * getTransformedNofe(egamma2, b)
0192                      * getSigma(egamma1) * getSigma(egamma2) *_probOfBreakup[k+1]*b;
0193                 bint = bint + 0.5*(f1+f2)*(b-bold);
0194             }
0195             bint = 2.0*starlightConstants::pi*bint;
0196             _fnDouble[i][j] = bint * (_eGamma[i+1] - _eGamma[i]) * (_eGamma[j+1] - _eGamma[j]);
0197             egamma2 = egamma2 * eg_inc;
0198         }
0199         egamma1 = egamma1 * eg_inc;
0200     }
0201
0202     for (int i = 0; i < _nK; i++)
0203     {
0204         _fnDoubleInt[i] = 0.0;
0205         for (int j = 0; j < _nK; j++)
0206         {
0207             _fnDoubleInt[i] = _fnDoubleInt[i] + _fnDouble[i][j];
0208         }
0209     }
0210
0211     _fnDoubleIntCumulative[0] = 0.0;
0212     for (int i = 1; i < _nK+1; i++)
0213     {
0214         _fnDoubleIntCumulative[i] = _fnDoubleIntCumulative[i-1] + _fnDoubleInt[i-1];
0215     }
0216
0217     fnorm = 1.0/_fnDoubleIntCumulative[_nK];
0218     for (int i = 0; i < _nK+1; i++)
0219     {
0220         _fnDoubleIntCumulative[i] = fnorm * _fnDoubleIntCumulative[i];
0221     }
0222
0223     return 0;
0224 }
0225
0226 double spectrum::drawKsingle()
0227 {
0228     double xtest = 0;
0229     int itest = 0;
0230     double egamma = 0.0;
0231
0232     xtest = _randy.Rndom();
0233     while (xtest > _fnSingleCumulative[itest])
0234     {
0235         itest++;
0236     }
0237     itest = itest - 1;
0238
0239     if (itest >= _nK || itest < 0)
0240     {
0241         std::cerr << "ERROR: itest: " << itest << std::endl;
0242
0243     }
0244
0245     double delta_f = xtest - _fnSingleCumulative[itest];
0246     if (delta_f <= 0.0)
0247     {
0248         std::cout << "WARNING: delta_f: " << delta_f << std::endl;
0249         return -1;
0250     }
0251     double dE = _eGamma[itest+1] - _eGamma[itest];
0252     double dF = _fnSingleCumulative[itest+1] - _fnSingleCumulative[itest];
0253
0254     double delta_e = delta_f*dE/dF;
0255
0256     if (delta_e > (_eGamma[itest+1] - _eGamma[itest]))
0257     {
0258         std::cerr << "ERROR: delta_E: " << delta_e << std::endl;
0259     }
0260
0261     egamma = _eGamma[itest] + delta_e;
0262     return egamma;
0263 }
0264
0265 void spectrum::drawKdouble(float& egamma1, float& egamma2)
0266 {
0267     double xtest1 = 0.0;
0268     double xtest2 = 0.0;
0269     int itest1 = 0;
0270     int itest2 = 0;
0271
0272     xtest1 = _randy.Rndom();
0273
0274     while (xtest1 > _fnDoubleIntCumulative[itest1])
0275     {
0276         itest1++;
0277     }
0278     itest1 = itest1 - 1;
0279
0280     if (itest1 >= _nK || itest1 < 0)
0281     {
0282         std::cerr << "ERROR: itest1: " << itest1 << std::endl;
0283     }
0284     double delta_f = xtest1 - _fnDoubleIntCumulative[itest1];
0285
0286     if (delta_f <= 0.0)
0287     {
0288         std::cout << "WARNING: delta_f: " << delta_f << std::endl;
0289     }
0290
0291     double dE = _eGamma[itest1+1] - _eGamma[itest1];
0292     double dF = _fnDoubleIntCumulative[itest1+1] - _fnDoubleIntCumulative[itest1];
0293
0294     double delta_e = delta_f*dE/dF;
0295
0296     if (delta_e > (_eGamma[itest1+1] - _eGamma[itest1]))
0297     {
0298         std::cerr << "ERROR: delta_E: " << delta_e << std::endl;
0299     }
0300
0301     egamma1 = _eGamma[itest1] + delta_e;
0302
0303     // Second gamma
0304
0305     std::vector<double> fn_second_cumulative(_nK+1);
0306
0307     fn_second_cumulative[0] = 0.0;
0308     for(int i = 1; i < _nK+1; i++)
0309     {
0310        fn_second_cumulative[i] = fn_second_cumulative[i-1] + _fnDouble[itest1][i-1];
0311     }
0312
0313     double norm_factor = 1.0/fn_second_cumulative[_nK];
0314     for(int i = 0; i < _nK+1; i++)
0315     {
0316       fn_second_cumulative[i] = norm_factor*fn_second_cumulative[i];
0317     }
0318
0319     xtest2 = _randy.Rndom();
0320
0321     while (xtest2 > fn_second_cumulative[itest2])
0322     {
0323         itest2++;
0324     }
0325     itest2 = itest2 - 1;
0326
0327     if (itest2 >= _nK || itest2 < 0)
0328     {
0329         std::cerr << "ERROR: itest2: " << itest2 << std::endl;
0330     }
0331     delta_f = xtest2 - fn_second_cumulative[itest2];
0332
0333     if (delta_f <= 0.0)
0334     {
0335         std::cout << "WARNING: delta_f: " << delta_f << std::endl;
0336     }
0337
0338     dE = _eGamma[itest2+1] - _eGamma[itest2];
0339     dF = fn_second_cumulative[itest2+1] - fn_second_cumulative[itest2];
0340
0341     delta_e = delta_f*dE/dF;
0342
0343     if (delta_e > (_eGamma[itest2+1] - _eGamma[itest2]))
0344     {
0345         std::cerr << "ERROR: delta_E: " << delta_e << std::endl;
0346     }
0347
0348     egamma2 = _eGamma[itest2] + delta_e;
0349
0350     return;
0351 }
0352
0353
0354 bool spectrum::generateBreakupProbabilities()
0355 {
0356
0357     int nbbins = _nBbins;
0358
0359     double b_min = _bMin;
0360     double binc = exp((log(_bMax/_bMin))/(double)_nBbins);
0361
0362     double b = b_min;
0363
0364     _probOfBreakup.resize(nbbins);
0365
0366     for (int i = 0; i < nbbins; i++)
0367     {
0368         double bimp = b;
0369         double rhad = 0;
0370         rhad = _beamBeamSystem->probabilityOfBreakup(bimp);
0371         _probOfBreakup[i] = exp(-rhad);
0372         b = b*binc;
0373     }
0374     return true;
0375 }
0376
0377 double spectrum::getFnSingle(double egamma) const
0378 {
0379     double eginc = exp(log(_eGammaMax/_eGammaMin)/static_cast<double>(_nK));
0380     int i1 = log(egamma/_eGammaMin)/log(eginc);
0381     int i2 = i1 + 1;
0382     double fnSingle = 0.0;
0383
0384     if (i1 < 0 || i2 > _nK)
0385     {
0386         std::cout << "I1, I2 out of bounds. Egamma = " << egamma << std::endl;
0387         std::cout << "I1, I2 = " << i1 << ", " << i2 << std::endl;
0388         fnSingle = 0.0;
0389     }
0390     else
0391     {
0392         double dE = _eGamma[i2] - _eGamma[i1];
0393         double eFrac = (egamma - _eGamma[i1])/dE;
0394
0395         if (eFrac < 0.0 || eFrac > 1.0)
0396         {
0397             std::cout << "WARNING: Efrac = " << eFrac << std::endl;
0398         }
0399         fnSingle = (1.0 - eFrac)*_fnSingle[i1] + eFrac*_fnSingle[i2];
0400     }
0401     return fnSingle;
0402 }
0403
0404 double spectrum::getFnDouble(double egamma1, double egamma2) const
0405 {
0406     double eginc = exp(log(_eGammaMax/_eGammaMin)/static_cast<double>(_nK));
0407     int i1 = log(egamma1/_eGammaMin)/log(eginc);
0408     int i2 = i1 + 1;
0409     double fnDouble = 0.0;
0410
0411     if (i1 < 0 || i2 > _nK)
0412     {
0413         std::cout << "I1, I2 out of bounds. Egamma1 = " << egamma1 << std::endl;
0414         std::cout << "I1, I2 = " << i1 << ", " << i2 << std::endl;
0415         fnDouble = 0.0;
0416         return fnDouble;
0417     }
0418
0419     int j1 = log(egamma2/_eGammaMin)/log(eginc);
0420     int j2 = j1 + 1;
0421
0422     if (j1 < 0 || j2 > _nK)
0423     {
0424         std::cout << "J1, J2 out of bounds. Egamma2 = " << egamma2 << std::endl;
0425         std::cout << "J1, J2 = " << j1 << ", " << j2 << std::endl;
0426         fnDouble = 0.0;
0427         return fnDouble;
0428     }
0429
0430     double dE1 = _eGamma[i2] - _eGamma[i1];
0431     double eFrac1 = (egamma1 - _eGamma[i1])/dE1;
0432
0433     if (eFrac1 < 0.0 || eFrac1 > 1.0)
0434     {
0435         std::cout << "WARNING: Efrac1 = " << eFrac1 << std::endl;
0436     }
0437
0438     double ptemp1 = (1.0 - eFrac1)*_fnDouble[i1][j1] + eFrac1*_fnDouble[i2][j1];
0439     double ptemp2 = (1.0 - eFrac1)*_fnDouble[i1][j2] + eFrac1*_fnDouble[i2][j2];
0440
0441     double dE2 = _eGamma[j2] - _eGamma[j1];
0442     double eFrac2 = (egamma2 - _eGamma[j1])/dE2;
0443
0444     if (eFrac2 < 0.0 || eFrac2 > 1.0)
0445     {
0446         std::cout << "WARNING: Efrac2 = " << eFrac2 << std::endl;
0447     }
0448
0449     fnDouble = (1.0 - eFrac2)*ptemp1 + eFrac2*ptemp2;
0450
0451     return fnDouble;
0452
0453 }
0454
0455 double spectrum::getTransformedNofe(double egamma, double b)
0456 {
0457    double factor = 1.0/(2.0*_beamBeamSystem->beamLorentzGamma());
0458    double res = factor * _beamBeamSystem->beam1().photonDensity(b, egamma*factor);
0459
0460    return res;
0461 }
0462
0463
0464
0465