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

   1
   2 /*
   3     <one line to give the program's name and a brief idea of what it does.>
   4     Copyright (C) <year>  <name of author>
   5
   6 p    This program is free software: you can redistribute it and/or modify
   7     it under the terms of the GNU General Public License as published by
   8     the Free Software Foundation, either version 3 of the License, or
   9     (at your option) any later version.
  10
  11     This program is distributed in the hope that it will be useful,
  12     but WITHOUT ANY WARRANTY; without even the implied warranty of
  13     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14     GNU General Public License for more details.
  15
  16     You should have received a copy of the GNU General Public License
  17     along with this program.  If not, see <http://www.gnu.org/licenses/>.
  18
  19 */
  20
  21 #include "spectrum.h"
  22 #include <cmath>
  23 #include "beambeamsystem.h"
  24 #include <randomgenerator.h>
  25 #include <iostream>
  26
  27 spectrum::spectrum(beamBeamSystem *bbs) :
  28          _bMin(5.0)
  29         ,_bMax(128000.0)
  30         ,_nBbins(6400)
  31         ,_probOfBreakup(_nBbins)
  32         ,_beamBeamSystem(bbs)
  33         ,_nK(10000)
  34         ,_fnSingle(_nK)
  35         ,_fnDouble(_nK)
  36         ,_fnSingleCumulative(_nK+1)
  37         ,_fnDoubleCumulative(_nK+1)
  38         ,_fnDoubleInt(_nK)
  39         ,_fnDoubleIntCumulative(_nK+1)
  40         ,_eGamma(_nK+1)
  41         ,_eGammaMin(6.0)
  42         ,_eGammaMax(600000.0)
  43         ,_zTarget(82)
  44         ,_aTarget(278)
  45         ,_hadBreakProbCalculated(false)
  46 {
  47     _eGamma.resize(_nK+1);
  48     _probOfBreakup.resize(_nBbins);
  49 }
  50
  51 int spectrum::generateKsingle()
  52 {
  53
  54     _fnSingle.resize(_nK);
  55     _fnSingleCumulative.resize(_nK+1);
  56
  57     double eg_inc = exp(log(_eGammaMax/_eGammaMin)/(double)_nK);
  58
  59     double egamma =  _eGammaMin;
  60     for (int i = 0; i < _nK+1; i++)
  61     {
  62         _eGamma[i] = egamma;
  63         egamma = egamma * eg_inc;
  64     }
  65     egamma = _eGammaMin;
  66
  67     double fnorm = 0;
  68
  69
  70     if (_hadBreakProbCalculated == false)
  71     {
  72         _hadBreakProbCalculated = generateBreakupProbabilities();
  73     }
  74     double binc = exp((log(_bMax/_bMin))/(double)_nBbins);
  75
  76     for (int i = 0; i < _nK; i++)
  77     {
  78         double b = _bMin;
  79
  80         double bint = 0.0;
  81
  82         double f1 = 0;
  83         double f2 = 0;
  84
  85         for (int j = 0; j < _nBbins - 1; j++)
  86         {
  87             double bold = b;
  88             if (j == 0)
  89             {
  90                 //f1 = fBeamBeamSystem->getBeam1().nofe(egamma, b)*GetSigma(egamma)*fProbOfBreakup[j]*b;
  91                 f1 = getTransformedNofe(egamma, b)*getSigma(egamma)*_probOfBreakup[j]*b;
  92                 //std::cout << fProbOfBreakup[j] << std::endl;
  93             }
  94             else
  95             {
  96                 f1 = f2;
  97             }
  98             b = b*binc;
  99 //            f2 = fBeamBeamSystem->getBeam1().nofe(egamma, b)*GetSigma(egamma)*fProbOfBreakup[j+1]*b;;
 100             f2 = getTransformedNofe(egamma, b)*getSigma(egamma)*_probOfBreakup[j+1]*b;;
 101             bint = bint + 0.5*(f1+f2)*(b-bold);
 102         }
 103         bint = 2.0*starlightConstants::pi*bint;
 104         if (i == 0)
 105         {
 106             fnorm = 1.0/bint;
 107         }
 108         _fnSingle[i] = bint*(_eGamma[i+1]-_eGamma[i]);
 109
 110         egamma = egamma*eg_inc;
 111     }
 112
 113     _fnSingleCumulative[0] = 0.00;
 114     for (int i = 0; i < _nK; i++)
 115     {
 116         _fnSingleCumulative[i+1] = _fnSingleCumulative[i]+_fnSingle[i];
 117     }
 118
 119     double fnormfactor = 1.00/_fnSingleCumulative[_nK];
 120     for (int i = 0; i < _nK; i++)
 121     {
 122         _fnSingleCumulative[i+1] = fnormfactor*_fnSingleCumulative[i+1];
 123     }
 124
 125     return 0;
 126
 127 }
 128
 129 int spectrum::generateKdouble()
 130 {
 131     //Quick fix for now TODO: Fix it!
 132     _nK = 100;
 133
 134     _fnDouble.resize(_nK);
 135     _fnDoubleInt.resize(_nK);
 136     _fnDoubleIntCumulative.resize(_nK+1);
 137     _fnDoubleCumulative.resize(_nK+1);
 138     for (int i = 0; i < _nK; i++)
 139     {
 140         _fnDouble[i].resize(_nK);
 141         _fnDoubleCumulative[i].resize(_nK+1);
 142     }
 143     _fnDoubleCumulative[_nK].resize(_nK+1);
 144
 145     double eg_inc = exp(log(_eGammaMax/_eGammaMin)/(double)_nK);
 146     double egamma1 =  _eGammaMin;
 147     double egamma2 =  _eGammaMin;
 148
 149     for (int i = 0; i < _nK+1; i++)
 150     {
 151         _eGamma[i] = egamma1;
 152         egamma1 = egamma1 * eg_inc;
 153     }
 154     egamma1 = _eGammaMin;
 155
 156     double fnorm = 0;
 157
 158     if (_hadBreakProbCalculated == false)
 159     {
 160         _hadBreakProbCalculated = generateBreakupProbabilities();
 161     }
 162
 163     double binc = exp((log(_bMax/_bMin))/(double)_nBbins);
 164
 165     int nbbins = _nBbins;
 166
 167     //double b_min = _bMin;
 168     //double b_max = _bMax;
 169
 170     for (int i = 0; i < _nK; i++)
 171     {
 172
 173         egamma2 = _eGammaMin;
 174         //double sum_over_k = 0.0;
 175
 176         for (int j = 0; j < _nK; j++)
 177         {
 178             double bint = 0.0;
 179             double b = _bMin;
 180             double f1 = 0;
 181             double f2 = 0;
 182
 183             for (int k = 0; k < nbbins - 1; k++)
 184             {
 185                 double bold = b;
 186
 187                 if (k == 0)
 188                 {
 189                   // f1 = fBeamBeamSystem->getBeam1().nofe(egamma1, b) * fBeamBeamSystem->getBeam2().nofe(egamma2, b)
 190                   //       * GetSigma(egamma1) * GetSigma(egamma2) *fProbOfBreakup[j]*b;
 191                   f1 = getTransformedNofe(egamma1, b) * getTransformedNofe(egamma2, b)
 192                          * getSigma(egamma1) * getSigma(egamma2) *_probOfBreakup[k]*b;                }
 193                 else
 194                 {
 195                     f1 = f2;
 196                 }
 197                 b = b*binc;
 198                 // f2 = fBeamBeamSystem->getBeam1().nofe(egamma1, b) * fBeamBeamSystem->getBeam2().nofe(egamma2, b)
 199                 //     * GetSigma(egamma1) * GetSigma(egamma2) *fProbOfBreakup[j+1]*b;
 200                 f2 = getTransformedNofe(egamma1, b) * getTransformedNofe(egamma2, b)
 201                      * getSigma(egamma1) * getSigma(egamma2) *_probOfBreakup[k+1]*b;
 202                 bint = bint + 0.5*(f1+f2)*(b-bold);
 203             }
 204             bint = 2.0*starlightConstants::pi*bint;
 205             _fnDouble[i][j] = bint * (_eGamma[i+1] - _eGamma[i]) * (_eGamma[j+1] - _eGamma[j]);
 206             egamma2 = egamma2 * eg_inc;
 207         }
 208         egamma1 = egamma1 * eg_inc;
 209     }
 210
 211     for (int i = 0; i < _nK; i++)
 212     {
 213         _fnDoubleInt[i] = 0.0;
 214         for (int j = 0; j < _nK; j++)
 215         {
 216             _fnDoubleInt[i] = _fnDoubleInt[i] + _fnDouble[i][j];
 217         }
 218     }
 219
 220     _fnDoubleIntCumulative[0] = 0.0;
 221     for (int i = 1; i < _nK+1; i++)
 222     {
 223         _fnDoubleIntCumulative[i] = _fnDoubleIntCumulative[i-1] + _fnDoubleInt[i-1];
 224     }
 225
 226     fnorm = 1.0/_fnDoubleIntCumulative[_nK];
 227     for (int i = 0; i < _nK+1; i++)
 228     {
 229         _fnDoubleIntCumulative[i] = fnorm * _fnDoubleIntCumulative[i];
 230     }
 231
 232     return 0;
 233 }
 234
 235 double spectrum::drawKsingle()
 236 {
 237     double xtest = 0;
 238     int itest = 0;
 239     double egamma = 0.0;
 240
 241     xtest = randyInstance.Rndom();
 242     while (xtest > _fnSingleCumulative[itest])
 243     {
 244         itest++;
 245     }
 246     itest = itest - 1;
 247
 248     if (itest >= _nK || itest < 0)
 249     {
 250         std::cerr << "ERROR: itest: " << itest << std::endl;
 251
 252     }
 253
 254     double delta_f = xtest - _fnSingleCumulative[itest];
 255     if (delta_f <= 0.0)
 256     {
 257         std::cout << "WARNING: delta_f: " << delta_f << std::endl;
 258         return -1;
 259     }
 260     double dE = _eGamma[itest+1] - _eGamma[itest];
 261     double dF = _fnSingleCumulative[itest+1] - _fnSingleCumulative[itest];
 262
 263     double delta_e = delta_f*dE/dF;
 264
 265     if (delta_e > (_eGamma[itest+1] - _eGamma[itest]))
 266     {
 267         std::cerr << "ERROR: delta_E: " << delta_e << std::endl;
 268     }
 269
 270     egamma = _eGamma[itest] + delta_e;
 271     return egamma;
 272 }
 273
 274 void spectrum::drawKdouble(float& egamma1, float& egamma2)
 275 {
 276     double xtest1 = 0.0;
 277     double xtest2 = 0.0;
 278     int itest1 = 0;
 279     int itest2 = 0;
 280
 281     xtest1 = randyInstance.Rndom();
 282
 283     while (xtest1 > _fnDoubleIntCumulative[itest1])
 284     {
 285         itest1++;
 286     }
 287     itest1 = itest1 - 1;
 288
 289     if (itest1 >= _nK || itest1 < 0)
 290     {
 291         std::cerr << "ERROR: itest1: " << itest1 << std::endl;
 292     }
 293     double delta_f = xtest1 - _fnDoubleIntCumulative[itest1];
 294
 295     if (delta_f <= 0.0)
 296     {
 297         std::cout << "WARNING: delta_f: " << delta_f << std::endl;
 298     }
 299
 300     double dE = _eGamma[itest1+1] - _eGamma[itest1];
 301     double dF = _fnDoubleIntCumulative[itest1+1] - _fnDoubleIntCumulative[itest1];
 302
 303     double delta_e = delta_f*dE/dF;
 304
 305     if (delta_e > (_eGamma[itest1+1] - _eGamma[itest1]))
 306     {
 307         std::cerr << "ERROR: delta_E: " << delta_e << std::endl;
 308     }
 309
 310     egamma1 = _eGamma[itest1] + delta_e;
 311
 312     // Second gamma
 313
 314     //    double reldw = delta_e/(fEGamma[itest1+1] - fEGamma[itest1]);
 315     //    std::vector<double> fn_second(fNK);
 316     std::vector<double> fn_second_cumulative(_nK+1);
 317
 318     //    for(int i = 0; i < fNK; i++)
 319     //    {
 320     //       fn_second[i] = fFnDouble[itest1][i] + (fFnDouble[itest1+1][i] - fFnDouble[itest1][i])*reldw;
 321     //    }
 322
 323     fn_second_cumulative[0] = 0.0;
 324     for(int i = 1; i < _nK+1; i++)
 325     {
 326        //       fn_second_cumulative[i] = fn_second_cumulative[i-1] + fn_second[i-1]; //TODO:check indexing
 327        fn_second_cumulative[i] = fn_second_cumulative[i-1] + _fnDouble[itest1][i-1];
 328     }
 329
 330     double norm_factor = 1.0/fn_second_cumulative[_nK];
 331     for(int i = 0; i < _nK+1; i++)
 332     {
 333       fn_second_cumulative[i] = norm_factor*fn_second_cumulative[i];
 334     }
 335
 336     xtest2 = randyInstance.Rndom();
 337
 338     while (xtest2 > fn_second_cumulative[itest2])
 339     {
 340         itest2++;
 341     }
 342     itest2 = itest2 - 1;
 343
 344     if (itest2 >= _nK || itest2 < 0)
 345     {
 346         std::cerr << "ERROR: itest2: " << itest2 << std::endl;
 347     }
 348     delta_f = xtest2 - fn_second_cumulative[itest2];
 349
 350     if (delta_f <= 0.0)
 351     {
 352         std::cout << "WARNING: delta_f: " << delta_f << std::endl;
 353     }
 354
 355     dE = _eGamma[itest2+1] - _eGamma[itest2];
 356     dF = fn_second_cumulative[itest2+1] - fn_second_cumulative[itest2];
 357
 358     delta_e = delta_f*dE/dF;
 359
 360     if (delta_e > (_eGamma[itest2+1] - _eGamma[itest2]))
 361     {
 362         std::cerr << "ERROR: delta_E: " << delta_e << std::endl;
 363     }
 364
 365     egamma2 = _eGamma[itest2] + delta_e;
 366
 367     return;
 368 }
 369
 370
 371 bool spectrum::generateBreakupProbabilities()
 372 {
 373
 374     int nbbins = _nBbins;
 375
 376     double b_min = _bMin;
 377     //double b_max = _bMax;
 378
 379     //    double binc = (log(b_max/b_min))/(double)nbbins;
 380     double binc = exp((log(_bMax/_bMin))/(double)_nBbins);
 381
 382     double b = b_min;
 383
 384     _probOfBreakup.resize(nbbins);
 385
 386     for (int i = 0; i < nbbins; i++)
 387     {
 388         double bimp = b;
 389         double rhad = 0;
 390         rhad = _beamBeamSystem->probabilityOfBreakup(bimp);
 391         _probOfBreakup[i] = exp(-rhad);
 392         b = b*binc;
 393     }
 394     return true;
 395 }
 396
 397 double spectrum::getFnSingle(double egamma) const
 398 {
 399     double eginc = exp(log(_eGammaMax/_eGammaMin)/static_cast<double>(_nK));
 400     int i1 = log(egamma/_eGammaMin)/log(eginc);
 401     int i2 = i1 + 1;
 402     double fnSingle = 0.0;
 403
 404     if (i1 < 0 || i2 > _nK)
 405     {
 406         std::cout << "I1, I2 out of bounds. Egamma = " << egamma << std::endl;
 407         std::cout << "I1, I2 = " << i1 << ", " << i2 << std::endl;
 408         fnSingle = 0.0;
 409     }
 410     else
 411     {
 412         double dE = _eGamma[i2] - _eGamma[i1];
 413         double eFrac = (egamma - _eGamma[i1])/dE;
 414
 415         if (eFrac < 0.0 || eFrac > 1.0)
 416         {
 417             std::cout << "WARNING: Efrac = " << eFrac << std::endl;
 418         }
 419         fnSingle = (1.0 - eFrac)*_fnSingle[i1] + eFrac*_fnSingle[i2];
 420     }
 421     return fnSingle;
 422 }
 423
 424 double spectrum::getFnDouble(double egamma1, double egamma2) const
 425 {
 426     double eginc = exp(log(_eGammaMax/_eGammaMin)/static_cast<double>(_nK));
 427     int i1 = log(egamma1/_eGammaMin)/log(eginc);
 428     int i2 = i1 + 1;
 429     double fnDouble = 0.0;
 430
 431     if (i1 < 0 || i2 > _nK)
 432     {
 433         std::cout << "I1, I2 out of bounds. Egamma1 = " << egamma1 << std::endl;
 434         std::cout << "I1, I2 = " << i1 << ", " << i2 << std::endl;
 435         fnDouble = 0.0;
 436         return fnDouble;
 437     }
 438
 439     int j1 = log(egamma2/_eGammaMin)/log(eginc);
 440     int j2 = j1 + 1;
 441
 442     if (j1 < 0 || j2 > _nK)
 443     {
 444         std::cout << "J1, J2 out of bounds. Egamma2 = " << egamma2 << std::endl;
 445         std::cout << "J1, J2 = " << j1 << ", " << j2 << std::endl;
 446         fnDouble = 0.0;
 447         return fnDouble;
 448     }
 449
 450     double dE1 = _eGamma[i2] - _eGamma[i1];
 451     double eFrac1 = (egamma1 - _eGamma[i1])/dE1;
 452
 453     if (eFrac1 < 0.0 || eFrac1 > 1.0)
 454     {
 455         std::cout << "WARNING: Efrac1 = " << eFrac1 << std::endl;
 456     }
 457
 458     double ptemp1 = (1.0 - eFrac1)*_fnDouble[i1][j1] + eFrac1*_fnDouble[i2][j1];
 459     double ptemp2 = (1.0 - eFrac1)*_fnDouble[i1][j2] + eFrac1*_fnDouble[i2][j2];
 460
 461     double dE2 = _eGamma[j2] - _eGamma[j1];
 462     double eFrac2 = (egamma2 - _eGamma[j1])/dE2;
 463
 464     if (eFrac2 < 0.0 || eFrac2 > 1.0)
 465     {
 466         std::cout << "WARNING: Efrac2 = " << eFrac2 << std::endl;
 467     }
 468
 469     fnDouble = (1.0 - eFrac2)*ptemp1 + eFrac2*ptemp2;
 470
 471     return fnDouble;
 472
 473 }
 474
 475 double spectrum::getTransformedNofe(double egamma, double b)
 476 {
 477    double factor = 1.0/(2.0*_beamBeamSystem->beamLorentzGamma());
 478    double res = factor * _beamBeamSystem->beam1().photonFlux(b, egamma*factor);
 479
 480    return res;
 481 }
 482
 483
 484
 485