STARLIGHT/starlight/src/.svn/text-base/twophotonluminosity.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 //    Added incoherent factor to luminosity table output--Joey
  29 //
  30 //
  31 ///////////////////////////////////////////////////////////////////////////
  32
  33
  34 #include <iostream>
  35 #include <fstream>
  36 #include <cmath>
  37
  38 #include "inputParameters.h"
  39 #include "beambeamsystem.h"
  40 #include "beam.h"
  41 #include "starlightconstants.h"
  42 #include "nucleus.h"
  43 #include "bessel.h"
  44 #include "twophotonluminosity.h"
  45 #include <pthread.h>
  46
  47 using namespace std;
  48 using namespace starlightConstants;
  49
  50
  51
  52 //______________________________________________________________________________
  53 twoPhotonLuminosity::twoPhotonLuminosity(beam beam_1,beam beam_2):
  54 beamBeamSystem(beam_1,beam_2)
  55 ,_gamma(inputParametersInstance.beamLorentzGamma())
  56 ,_nWbins(inputParametersInstance.nmbWBins())
  57 ,_nYbins(inputParametersInstance.nmbRapidityBins())
  58 ,_wMin(inputParametersInstance.minW())
  59 ,_yMin(-inputParametersInstance.maxRapidity())
  60 ,_wMax(inputParametersInstance.maxW())
  61 ,_yMax(inputParametersInstance.maxRapidity())
  62 {
  63   //Lets check to see if we need to recalculate the luminosity tables
  64   twoPhotonDifferentialLuminosity();
  65 }
  66
  67
  68 //______________________________________________________________________________
  69 twoPhotonLuminosity::~twoPhotonLuminosity()
  70 { }
  71
  72
  73 //______________________________________________________________________________
  74 void twoPhotonLuminosity::twoPhotonDifferentialLuminosity()
  75 {
  76   ofstream wylumfile;
  77   wylumfile.precision(15);
  78   wylumfile.open("slight.txt");
  79   std::vector<double> w(_nWbins);
  80   std::vector<double> y(_nYbins);
  81   double xlum = 0.;
  82   double Normalize = 0.,OldNorm;
  83   double wmev = 0;
  84
  85   Normalize = 1./sqrt(1*(double)_nWbins*_nYbins); //if your grid is very fine, you'll want high accuracy-->small Normalize
  86   OldNorm   = Normalize;
  87
  88   //Writing out our input parameters+(w,y)grid+diff._lum.
  89   wylumfile << inputParametersInstance.parameterValueKey() << endl;
  90   wylumfile << beam1().Z() <<endl;
  91   wylumfile << beam1().A() <<endl;
  92   wylumfile << beam2().Z() <<endl;
  93   wylumfile << beam2().A() <<endl;
  94   wylumfile << inputParametersInstance.beamLorentzGamma() <<endl;
  95   wylumfile << inputParametersInstance.maxW() <<endl;
  96   wylumfile << inputParametersInstance.minW() <<endl;
  97   wylumfile << inputParametersInstance.nmbWBins() <<endl;
  98   wylumfile << inputParametersInstance.maxRapidity() <<endl;
  99   wylumfile << inputParametersInstance.nmbRapidityBins() <<endl;
 100   wylumfile << inputParametersInstance.productionMode() <<endl;
 101   wylumfile << inputParametersInstance.beamBreakupMode() <<endl;
 102   wylumfile << inputParametersInstance.interferenceEnabled() <<endl;
 103   wylumfile << inputParametersInstance.interferenceStrength() <<endl;
 104   wylumfile << inputParametersInstance.coherentProduction() <<endl;
 105   wylumfile << inputParametersInstance.incoherentFactor() <<endl;
 106   wylumfile << inputParametersInstance.deuteronSlopePar() <<endl;
 107   wylumfile << inputParametersInstance.maxPtInterference() <<endl;
 108   wylumfile << inputParametersInstance.nmbPtBinsInterference() <<endl;
 109   for (unsigned int i = 0; i < _nWbins; ++i) {
 110     w[i] = _wMin + (_wMax-_wMin)/_nWbins*i;
 111     wylumfile << w[i] <<endl;
 112   }
 113   for (unsigned int i = 0; i < _nYbins; ++i) {
 114     y[i] = -_yMax + 2.*_yMax*i/(_nYbins);
 115     wylumfile << y[i] <<endl;
 116   }
 117
 118   if(inputParametersInstance.xsecCalcMethod() == 0) {
 119
 120     for (unsigned int i = 0; i < _nWbins; ++i) {   //For each (w,y) pair, calculate the diff. _lum
 121       for (unsigned int j = 0; j < _nYbins; ++j) {
 122         wmev = w[i]*1000.;
 123         xlum = wmev * D2LDMDY(wmev,y[j],Normalize);   //Convert photon flux dN/dW to Lorentz invariant photon number WdN/dW
 124         if (j==0) OldNorm = Normalize;       //Save value of integral for each new W(i) and Y(i)
 125         wylumfile << xlum <<endl;
 126       }
 127       Normalize = OldNorm;
 128   }
 129
 130   }
 131   else if(inputParametersInstance.xsecCalcMethod() == 1) {
 132
 133         const int nthreads = inputParametersInstance.nThreads();
 134         pthread_t threads[nthreads];
 135         difflumiargs args[nthreads];
 136
 137         for(int t = 0; t < nthreads; t++)
 138         {
 139             args[t].self = this;
 140         }
 141         for (unsigned int i = 1; i <= _nWbins; ++i) {   //For each (w,y) pair, calculate the diff. _lum
 142             printf("Calculating cross section: %2.0f %% \r", float(i)/float(_nWbins)*100);
 143             fflush(stdout);
 144             unsigned int r = 1;
 145             for(unsigned int j = 0; j < _nYbins/nthreads; ++j)
 146             {
 147
 148                 for(int t = 0; t < nthreads; t++)
 149                 {
 150                     args[t].m = w[i];
 151                     args[t].y = y[r];
 152
 153                     pthread_create(&threads[t], NULL, &twoPhotonLuminosity::D2LDMDY_Threaded, &args[t]);
 154                     r++;
 155                 }
 156                 for(int t = 0; t < nthreads; t++)
 157                 {
 158                     pthread_join(threads[t], NULL);
 159                     xlum = w[i] * args[t].res;
 160                     wylumfile << xlum <<endl;
 161                 }
 162             }
 163             for(unsigned int t = 0; t < _nYbins%nthreads; t++)
 164             {
 165                 args[t].m = w[i];
 166                 args[t].y = y[r];
 167
 168                 pthread_create(&threads[t], NULL, &twoPhotonLuminosity::D2LDMDY_Threaded, &args[t]);
 169                 r++;
 170             }
 171             for(unsigned int t = 0; t < _nYbins%nthreads; t++)
 172             {
 173                 pthread_join(threads[t], NULL);
 174                 xlum = w[i] * args[t].res;
 175                 wylumfile << xlum <<endl;
 176             }
 177         }
 178     }
 179
 180     wylumfile << inputParametersInstance.parameterValueKey() << endl;
 181     wylumfile.close();
 182     return;
 183 }
 184
 185 //______________________________________________________________________________
 186 double twoPhotonLuminosity::D2LDMDY(double M,double Y,double &Normalize)
 187 {
 188   // double differential luminosity
 189
 190   double D2LDMDYx = 0.;
 191
 192   _W1    =  M/2.0*exp(Y);
 193   _W2    =  M/2.0*exp(-Y);
 194   int Zin=beam1().Z();
 195   D2LDMDYx = 2.0/M*Zin*Zin*Zin*Zin*(starlightConstants::alpha*starlightConstants::alpha)*integral(Normalize);  //treats it as a symmetric collision
 196   Normalize = D2LDMDYx*M/(2.0*beam1().Z()*beam1().Z()*
 197                           beam1().Z()*beam1().Z()*
 198                           starlightConstants::alpha*starlightConstants::alpha);
 199   //Normalization also treats it as symmetric
 200   return D2LDMDYx;
 201 }
 202
 203
 204
 205 //______________________________________________________________________________
 206 double twoPhotonLuminosity::D2LDMDY(double M, double Y) const
 207 {
 208   // double differential luminosity
 209
 210   double D2LDMDYx = 0.;
 211   double w1    =  M/2.0*exp(Y);
 212   double w2    =  M/2.0*exp(-Y);
 213
 214   //int Z1=beam1().Z();
 215   //int Z2=beam2().Z();
 216
 217   double r_nuc1 = beam1().nuclearRadius();
 218   double r_nuc2 = beam2().nuclearRadius();
 219
 220   double b1min = r_nuc1;
 221   double b2min = r_nuc2;
 222
 223   double b1max = max(5.*_gamma*hbarc/w1,5*r_nuc1);
 224   double b2max = max(5.*_gamma*hbarc/w2,5*r_nuc2);
 225
 226   const int nbins_b1 = 120;
 227   const int nbins_b2 = 120;
 228
 229   double log_delta_b1 = (log(b1max)-log(b1min))/nbins_b1;
 230   double log_delta_b2 = (log(b2max)-log(b2min))/nbins_b2;
 231   double sum = 0;
 232   for(int i = 0; i < nbins_b1; ++i)
 233   {
 234       // Sum from nested integral
 235       double sum_b2 = 0;
 236       double b1_low = b1min*exp(i*log_delta_b1);
 237       double b1_high = b1min*exp((i+1)*log_delta_b1);
 238       double b1_cent = (b1_high+b1_low)/2.;
 239       for(int j = 0; j < nbins_b2; ++j)
 240       {
 241         // Sum from nested
 242         double sum_phi = 0;
 243         double b2_low = b2min*exp(j*log_delta_b2);
 244         double b2_high = b2min*exp((j+1)*log_delta_b2);
 245         double b2_cent = (b2_high+b2_low)/2.;
 246
 247         // Gaussian integration n = 10
 248         // Since cos is symmetric around 0 we only need 5 of the
 249         // points in the gaussian integration.
 250         const int ngi = 5;
 251         double weights[ngi] =
 252         {
 253           0.2955242247147529,
 254           0.2692667193099963,
 255           0.2190863625159820,
 256           0.1494513491505806,
 257           0.0666713443086881,
 258         };
 259
 260         double abscissas[ngi] =
 261         {
 262           -0.1488743389816312,
 263           -0.4333953941292472,
 264           -0.6794095682990244,
 265           -0.8650633666889845,
 266           -0.9739065285171717,
 267         };
 268
 269         for(int k = 0; k < ngi; ++k)
 270         {
 271             double b_rel = sqrt(b1_cent*b1_cent+b2_cent*b2_cent + 2.*b1_cent*b2_cent*cos(pi*(abscissas[k]+1)));
 272
 273             sum_phi += weights[k] * probabilityOfBreakup(b_rel)*2;
 274         }
 275         sum_b2 += beam2().photonFlux(b2_cent,w2)*pi*sum_phi*b2_cent*(b2_high-b2_low);
 276       }
 277
 278       sum += beam1().photonFlux(b1_cent, w1)*sum_b2*b1_cent*(b1_high-b1_low);
 279
 280   }
 281   D2LDMDYx = 2.*pi*M/2.*sum;
 282   return D2LDMDYx;
 283 }
 284
 285
 286 void * twoPhotonLuminosity::D2LDMDY_Threaded(void * a)
 287 {
 288   difflumiargs *args = (difflumiargs*)a;
 289   double M = args->m;
 290   double Y = args->y;
 291   args->res = args->self->D2LDMDY(M, Y);
 292
 293   return NULL;
 294 }
 295
 296   /*
 297   D2LDMDYx = 2.0/M*Zin*Zin*Zin*Zin*(starlightConstants::alpha*starlightConstants::alpha)*integral(Normalize);  //treats it as a symmetric collision
 298   Normalize = D2LDMDYx*M/(2.0*beam1().Z()*beam1().Z()*
 299                           beam1().Z()*beam1().Z()*
 300                           starlightConstants::alpha*starlightConstants::alpha);
 301   //Normalization also treats it as symmetric
 302   return D2LDMDYx;
 303
 304   */
 305
 306
 307 //______________________________________________________________________________
 308 double twoPhotonLuminosity::integral(double Normalize)
 309 {
 310   int NIter = 0;
 311   int NIterMin = 0;
 312   double EPS = 0.;
 313   double RM = 0.;
 314   double u1 = 0.;
 315   double u2 = 0.;
 316   double B1 = 0.;
 317   double B2 = 0.;
 318   double Integrala = 0.;
 319   double totsummary = 0.;
 320   double NEval      = 0.;
 321   double Lower[3];
 322   double Upper[3];
 323   double WK[500000]; //used to be [1000000]
 324   double Result, Summary, ResErr, NFNEVL;
 325
 326   EPS = .01*Normalize;   //This is EPS for integration, 1% of previous integral value.
 327   // Change this to the Woods-Saxon radius to be consistent with the older calculations (JN 230710)
 328   RM  = beam1().nuclearRadius()/starlightConstants::hbarcmev;  //Assumes symmetry?
 329   // RM  = beam1().woodSaxonRadius()/starlightConstants::hbarcmev;
 330
 331   NIter = 10000 + (int)1000000*(int)Normalize; //if integral value is very small, we don't do too many intertions to get precision down to 1%
 332   NIterMin = 600;
 333   u1 = 9.*_gamma/_W1; //upper boundary in B1
 334   u2 = 9.*_gamma/_W2; //upper boundary in B2
 335   B1 = .4*_gamma/_W1; //intermediate boundary in B1
 336   B2 = .4*_gamma/_W2; //intermediate boundary in B2
 337   //The trick is that u1,2 and b1,2 could be less than RM-the lower integration boundary, thus integration area splits into 4,2 or 1 pieces
 338
 339   if (u1 < RM){
 340     Integrala = 0;
 341     totsummary = 0;
 342     NEval      = 0;
 343   }
 344   else if (B1 > RM){
 345     if (u2 < RM){
 346       Integrala = 0;
 347       totsummary = 0;
 348       NEval      = 0;
 349     }
 350     else if (B2 > RM){            //integral has 4 parts
 351       Integrala = 0;
 352       totsummary = 40000;
 353       NEval      = 0;
 354       Lower[0]   = RM;       //1
 355       Lower[1]   = RM;       //2
 356       Lower[2]   = 0.;       //3
 357       Upper[2]   = 2.*starlightConstants::pi;    //3
 358       Upper[0]   = B1;       //1
 359       Upper[1]   = B2;       //2
 360       radmul(3,Lower,Upper,NIterMin,NIter,EPS,WK,NIter,Result,ResErr,NFNEVL,Summary);
 361       Integrala   = Integrala + Result;
 362       totsummary = totsummary + 1000*Summary;
 363       NEval      = NEval + NFNEVL;
 364       Upper[0]   = u1;       //1
 365       Upper[1]   = B2;       //2
 366       Lower[0]   = B1;       //1
 367       Lower[1]   = RM;       //2
 368       radmul(3,Lower,Upper,NIterMin,NIter,EPS,WK,NIter,Result,ResErr,NFNEVL,Summary);
 369       Integrala   = Integrala + Result;
 370       totsummary = totsummary + 100*Summary;
 371       NEval      = NEval + NFNEVL;
 372       Upper[0]   = B1;       //1
 373       Upper[1]   = u2;       //2
 374       Lower[0]   = RM;       //1
 375       Lower[1]   = B2;       //2
 376       radmul(3,Lower,Upper,NIterMin,NIter,EPS,WK,NIter,Result,ResErr,NFNEVL,Summary);
 377       Integrala   = Integrala + Result;
 378       totsummary = totsummary + 100*Summary;
 379       NEval      = NEval + NFNEVL;
 380       Upper[0]   = u1;       //1
 381       Upper[1]   = u2;       //2
 382       Lower[0]   = B1;       //1
 383       Lower[1]   = B2;       //2
 384       radmul(3,Lower,Upper,NIterMin,NIter,EPS,WK,NIter,Result,ResErr,NFNEVL,Summary);
 385       Integrala   = Integrala + Result;
 386       totsummary = totsummary + Summary;
 387       NEval      = NEval + NFNEVL;
 388     }
 389     else {
 390       //integral has 2 parts, b2 integral has only 1 component
 391       Integrala   = 0;
 392       totsummary = 20000;
 393       NEval      = 0;
 394       Lower[0]   = RM;       //1
 395       Lower[1]   = RM;       //2
 396       Lower[2]   = 0.;       //3
 397       Upper[2]   = 2.*starlightConstants::pi;    //3
 398       Upper[0]   = B1;       //1
 399       Upper[1]   = u2;       //2
 400       radmul(3,Lower,Upper,NIterMin,NIter,EPS,WK,NIter,Result,ResErr,NFNEVL,Summary);
 401       Integrala   = Integrala + Result;
 402       totsummary = totsummary + 100*Summary;
 403       NEval      = NEval + NFNEVL;
 404       Upper[0]   = u1;       //1
 405       Lower[0]   = B1;       //1
 406       radmul(3,Lower,Upper,NIterMin,NIter,EPS,WK,NIter,Result,ResErr,NFNEVL,Summary);
 407       Integrala   = Integrala + Result;
 408       totsummary = totsummary + Summary;
 409       NEval      = NEval + NFNEVL;
 410     }
 411   }
 412   else{
 413     if (u2 < RM ){
 414       Integrala   = 0;
 415       totsummary = 0;
 416       NEval      = 0;
 417     }
 418     else if (B2 > RM){
 419       //integral has 2 parts, b1 integral has only 1 component
 420       Integrala   = 0;
 421       totsummary = 20000;
 422       NEval      = 0;
 423       Lower[0]   = RM;       //1
 424       Lower[1]   = RM;       //2
 425       Lower[2]   = 0.;       //2
 426       Upper[2]   = 2.*starlightConstants::pi;    //3
 427       Upper[0]   = u1;       //1
 428       Upper[1]   = B2;       //2
 429       radmul(3,Lower,Upper,NIterMin,NIter,EPS,WK,NIter,Result,ResErr,NFNEVL,Summary);
 430       Integrala   = Integrala + Result;
 431       totsummary = totsummary + 100*Summary;
 432       NEval      = NEval + NFNEVL;
 433       Upper[1]   = u2;       //2
 434       Lower[1]   = B2;       //2
 435       radmul(3,Lower,Upper,NIterMin,NIter,EPS,WK,NIter,Result,ResErr,NFNEVL,Summary);
 436       Integrala   = Integrala + Result;
 437       totsummary = totsummary + Summary;
 438       NEval      = NEval + NFNEVL;
 439     }
 440     else{                 //integral has 1 part
 441       Integrala   = 0;
 442       totsummary = 10000;
 443       NEval      = 0;
 444       Lower[0]   = RM;       //1
 445       Lower[1]   = RM;       //2
 446       Lower[2]   = 0.;       //3
 447       Upper[2]   = 2.*starlightConstants::pi;    //3
 448       Upper[0]   = u1;       //1
 449       Upper[1]   = u2;       //2
 450       radmul(3,Lower,Upper,NIterMin,NIter,EPS,WK,NIter,Result,ResErr,NFNEVL,Summary);
 451       Integrala   = Integrala + Result;
 452       totsummary = totsummary + Summary;
 453       NEval      = NEval + NFNEVL;
 454     }
 455   }
 456   Integrala = 2*starlightConstants::pi*Integrala;
 457   return Integrala;
 458 }
 459
 460 //______________________________________________________________________________
 461 double twoPhotonLuminosity::radmul(int N,double *A,double *B,int MINPTS,int MAXPTS,double EPS,double *WK,int IWK,double &RESULT,double &RELERR,double &NFNEVL,double &IFAIL)
 462 {
 463   double wn1[14] = {       -0.193872885230909911, -0.555606360818980835,
 464                            -0.876695625666819078, -1.15714067977442459,  -1.39694152314179743,
 465                            -1.59609815576893754,  -1.75461057765584494,  -1.87247878880251983,
 466                            -1.94970278920896201,  -1.98628257887517146,  -1.98221815780114818,
 467                            -1.93750952598689219,  -1.85215668343240347,  -1.72615963013768225};
 468
 469   double wn3[14] = {        0.0518213686937966768,    0.0314992633236803330,
 470                             0.0111771579535639891, -0.00914494741655235473,  -0.0294670527866686986,
 471                             -0.0497891581567850424, -0.0701112635269013768,   -0.0904333688970177241,
 472                             -0.110755474267134071,  -0.131077579637250419,    -0.151399685007366752,
 473                             -0.171721790377483099,  -0.192043895747599447,    -0.212366001117715794};
 474
 475     double wn5[14] = {        0.871183254585174982e-01,  0.435591627292587508e-01,
 476                               0.217795813646293754e-01, 0.108897906823146873e-01,  0.544489534115734364e-02,
 477                               0.272244767057867193e-02, 0.136122383528933596e-02,  0.680611917644667955e-03,
 478                               0.340305958822333977e-03, 0.170152979411166995e-03,  0.850764897055834977e-04,
 479                               0.425382448527917472e-04, 0.212691224263958736e-04,  0.106345612131979372e-04};
 480
 481     double wpn1[14] = {       -1.33196159122085045, -2.29218106995884763,
 482                               -3.11522633744855959, -3.80109739368998611, -4.34979423868312742,
 483                               -4.76131687242798352, -5.03566529492455417, -5.17283950617283939,
 484                               -5.17283950617283939, -5.03566529492455417, -4.76131687242798352,
 485                               -4.34979423868312742, -3.80109739368998611, -3.11522633744855959};
 486
 487     double wpn3[14] = {        0.0445816186556927292, -0.0240054869684499309,
 488                                -0.0925925925925925875, -0.161179698216735251,  -0.229766803840877915,
 489                                -0.298353909465020564,  -0.366941015089163228,  -0.435528120713305891,
 490                                -0.504115226337448555,  -0.572702331961591218,  -0.641289437585733882,
 491                                -0.709876543209876532,  -0.778463648834019195,  -0.847050754458161859};
 492
 493     RESULT = 0;
 494
 495     double ABSERR = 0.;
 496     double ctr[15], wth[15], wthl[15], z[15];
 497     double R1  = 1.;
 498     double HF  = R1/2.;
 499     double xl2 = 0.358568582800318073;
 500     double xl4 = 0.948683298050513796;
 501     double xl5 = 0.688247201611685289;
 502     double w2  = 980./6561.;
 503     double w4  = 200./19683.;
 504     double wp2 = 245./486.;
 505     double wp4 = 25./729.;
 506     int j1 =0;
 507     double SUM1, SUM2, SUM3, SUM4, SUM5, RGNCMP=0., RGNVAL, RGNERR, F2, F3, DIF, DIFMAX, IDVAXN =0.;
 508     IFAIL  = 3;
 509     if (N < 2 || N > 15) return 0;
 510     if (MINPTS > MAXPTS) return 0;
 511     int IFNCLS = 0;
 512     int IDVAX0 = 0;
 513     bool LDV = false;
 514     double TWONDM = pow(2.,(double)(N));
 515     int IRGNST = 2*N+3;
 516     int IRLCLS = (int)pow(2.,(N))+2*N*(N+1)+1;
 517     int ISBRGN = IRGNST;
 518     int ISBTMP, ISBTPP;
 519     int ISBRGS = IRGNST;
 520
 521     if ( MAXPTS < IRLCLS ) return 0;
 522     for ( int j = 0; j < N; j++ ){  //10
 523       ctr[j]=(B[j] + A[j])*HF;
 524       wth[j]=(B[j] - A[j])*HF;
 525     }
 526  L20:
 527     double RGNVOL = TWONDM; //20
 528     for ( int j = 0; j < N; j++ ){            //30
 529       RGNVOL = RGNVOL*wth[j];
 530       z[j] = ctr[j];
 531     }
 532
 533     SUM1 = integrand(N,z);
 534
 535     DIFMAX = 0;
 536     SUM2   = 0;
 537     SUM3   = 0;
 538     for ( int j = 0; j < N; j++ ) {   //40
 539       z[j]=ctr[j]-xl2*wth[j];
 540       F2=integrand(N,z);
 541
 542       z[j]=ctr[j]+xl2*wth[j];
 543       F2=F2+integrand(N,z);
 544       wthl[j]=xl4*wth[j];
 545
 546       z[j]=ctr[j]-wthl[j];
 547       F3=integrand(N,z);
 548
 549       z[j]=ctr[j]+wthl[j];
 550       F3=F3+integrand(N,z);
 551       SUM2=SUM2+F2;
 552       SUM3=SUM3+F3;
 553       DIF=fabs(7.*F2-F3-12.*SUM1);
 554       DIFMAX=max(DIF,DIFMAX);
 555
 556       if ( DIFMAX == DIF) IDVAXN = j+1;
 557       z[j]=ctr[j];
 558
 559     }
 560
 561     SUM4   = 0;
 562     for ( int j = 1; j < N; j++){ //70
 563
 564         j1=j-1;
 565         for ( int k = j; k < N; k++){  //60
 566           for ( int l = 0; l < 2; l++){      //50
 567             wthl[j1]=-wthl[j1];
 568             z[j1]=ctr[j1]+wthl[j1];
 569             for ( int m = 0; m < 2; m++){   //50
 570               wthl[k]=-wthl[k];
 571               z[k]=ctr[k]+wthl[k];
 572               SUM4=SUM4+integrand(N,z);
 573             }
 574           }
 575           z[k]=ctr[k];
 576         }
 577         z[j1]=ctr[j1];
 578     }
 579
 580     SUM5  = 0;
 581
 582     for ( int j = 0; j < N; j++){             //80
 583       wthl[j]=-xl5*wth[j];
 584       z[j]=ctr[j]+wthl[j];
 585     }
 586  L90:
 587     SUM5=SUM5+integrand(N,z);   //line 90
 588
 589     for (int j = 0; j < N; j++){ //100
 590       wthl[j]=-wthl[j];
 591       z[j]=ctr[j]+wthl[j];
 592       if ( wthl[j] > 0. ) goto L90;
 593     }
 594
 595     RGNCMP = RGNVOL*(wpn1[N-2]*SUM1+wp2*SUM2+wpn3[N-2]*SUM3+wp4*SUM4);
 596     RGNVAL = wn1[N-2]*SUM1+w2*SUM2+wn3[N-2]*SUM3+w4*SUM4+wn5[N-2]*SUM5;
 597
 598     RGNVAL = RGNVOL*RGNVAL;
 599     RGNERR = fabs(RGNVAL-RGNCMP);
 600     RESULT = RESULT+RGNVAL;
 601     ABSERR = ABSERR+RGNERR;
 602     IFNCLS = IFNCLS+IRLCLS;
 603
 604
 605     if (LDV){
 606     L110:
 607
 608       ISBTMP = 2*ISBRGN;
 609
 610       if ( ISBTMP > ISBRGS ) goto L160;
 611       if ( ISBTMP < ISBRGS ){
 612         ISBTPP = ISBTMP + IRGNST;
 613
 614         if ( WK[ISBTMP-1] < WK[ISBTPP-1] ) ISBTMP = ISBTPP;
 615       }
 616
 617       if ( RGNERR >= WK[ISBTMP-1] ) goto L160;
 618       for ( int k = 0; k < IRGNST; k++){
 619         WK[ISBRGN-k-1] = WK[ISBTMP-k-1];
 620                 }
 621       ISBRGN = ISBTMP;
 622       goto L110;
 623     }
 624  L140:
 625
 626     ISBTMP = (ISBRGN/(2*IRGNST))*IRGNST;
 627
 628     if ( ISBTMP >= IRGNST && RGNERR > WK[ISBTMP-1] ){
 629       for ( int k = 0; k < IRGNST; k++){
 630         WK[ISBRGN-k-1]=WK[ISBTMP-k-1];
 631       }
 632       ISBRGN = ISBTMP;
 633       goto L140;
 634     }
 635  L160:
 636
 637     WK[ISBRGN-1] = RGNERR;
 638     WK[ISBRGN-2] = RGNVAL;
 639     WK[ISBRGN-3] = IDVAXN;
 640
 641     for ( int j = 0; j < N; j++) {
 642
 643       ISBTMP = ISBRGN-2*j-4;
 644       WK[ISBTMP]=ctr[j];
 645       WK[ISBTMP-1]=wth[j];
 646     }
 647     if (LDV) {
 648       LDV = false;
 649       ctr[IDVAX0-1]=ctr[IDVAX0-1]+2*wth[IDVAX0-1];
 650       ISBRGS = ISBRGS + IRGNST;
 651       ISBRGN = ISBRGS;
 652       goto L20;
 653     }
 654
 655     RELERR=ABSERR/fabs(RESULT);
 656
 657
 658     if ( ISBRGS + IRGNST > IWK ) IFAIL = 2;
 659     if ( IFNCLS + 2*IRLCLS > MAXPTS ) IFAIL = 1;
 660     if ( RELERR < EPS && IFNCLS >= MINPTS ) IFAIL = 0;
 661
 662     if ( IFAIL == 3 ) {
 663       LDV = true;
 664       ISBRGN = IRGNST;
 665       ABSERR = ABSERR-WK[ISBRGN-1];
 666       RESULT = RESULT-WK[ISBRGN-2];
 667       IDVAX0 = (int)WK[ISBRGN-3];
 668
 669       for ( int j = 0; j < N; j++) {
 670         ISBTMP = ISBRGN-2*j-4;
 671         ctr[j] = WK[ISBTMP];
 672         wth[j] = WK[ISBTMP-1];
 673       }
 674
 675       wth[IDVAX0-1] = HF*wth[IDVAX0-1];
 676       ctr[IDVAX0-1] = ctr[IDVAX0-1]-wth[IDVAX0-1];
 677       goto L20;
 678     }
 679     NFNEVL=IFNCLS;
 680     return 1;
 681 }
 682
 683
 684 //______________________________________________________________________________
 685 double twoPhotonLuminosity::integrand(double ,  // N (unused)
 686                                       double X[])
 687 {
 688   double  b1 = X[0];      //1
 689   double  b2 = X[1];      //2
 690   double  theta = X[2];   //3
 691   //breakup effects distances in fermis, so convert to fermis(factor of hbarcmev)
 692   double  D = sqrt(b1*b1+b2*b2-2*b1*b2*cos(theta))*starlightConstants::hbarcmev;
 693   double  integrandx = Nphoton(_W1,_gamma,b1)*Nphoton(_W2,_gamma,b2)*b1*b2*probabilityOfBreakup(D);
 694   //why not just use gamma?
 695   //switching _input2photon.beamLorentzGamma()to gamma
 696   return integrandx;
 697 }
 698
 699
 700 //______________________________________________________________________________
 701 double twoPhotonLuminosity::Nphoton(double W,double gamma,double Rho)
 702 {
 703  double Nphoton1 =0.;
 704  double WGamma = W/gamma;
 705  double WGR    = 1.0*WGamma*Rho;
 706  //factor of Z^2*alpha is omitted
 707  double Wphib = WGamma*bessel::dbesk1(WGR);
 708
 709  Nphoton1 = 1.0/(starlightConstants::pi*starlightConstants::pi)*(Wphib*Wphib);
 710  return Nphoton1;
 711 }