estarlight/src/twophotonluminosity.cpp

0001 ///////////////////////////////////////////////////////////////////////////
0002 //
0003 //    Copyright 2010
0004 //
0005 //    This file is part of starlight.
0006 //
0007 //    starlight is free software: you can redistribute it and/or modify
0008 //    it under the terms of the GNU General Public License as published by
0009 //    the Free Software Foundation, either version 3 of the License, or
0010 //    (at your option) any later version.
0011 //
0012 //    starlight is distributed in the hope that it will be useful,
0013 //    but WITHOUT ANY WARRANTY; without even the implied warranty of
0014 //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0015 //    GNU General Public License for more details.
0016 //
0017 //    You should have received a copy of the GNU General Public License
0018 //    along with starlight. If not, see <http://www.gnu.org/licenses/>.
0019 //
0020 ///////////////////////////////////////////////////////////////////////////
0021 //
0022 // File and Version Information:
0023 // $Rev:: 271                         $: revision of last commit
0024 // $Author:: jnystrand                $: author of last commit
0025 // $Date:: 2016-07-08 17:28:57 +0100 #$: date of last commit
0026 //
0027 // Description:
0028 //    Added incoherent factor to luminosity table output--Joey
0029 //
0030 //
0031 ///////////////////////////////////////////////////////////////////////////
0032
0033
0034 #include <iostream>
0035 #include <fstream>
0036 #include <cmath>
0037
0038 #include "inputParameters.h"
0039 #include "beambeamsystem.h"
0040 #include "beam.h"
0041 #include "starlightconstants.h"
0042 #include "nucleus.h"
0043 #include "bessel.h"
0044 #include "twophotonluminosity.h"
0045
0046 using namespace std;
0047 using namespace starlightConstants;
0048
0049
0050
0051 //______________________________________________________________________________
0052 twoPhotonLuminosity::twoPhotonLuminosity(const inputParameters& inputParametersInstance, beam beam_1,beam beam_2):
0053 beamBeamSystem(inputParametersInstance, beam_1,beam_2)
0054 ,_gamma(inputParametersInstance.beamLorentzGamma())
0055 ,_nWbins(inputParametersInstance.nmbWBins())
0056 ,_nYbins(inputParametersInstance.nmbRapidityBins())
0057 ,_wMin(inputParametersInstance.minW())
0058 ,_yMin(-inputParametersInstance.maxRapidity())
0059 ,_wMax(inputParametersInstance.maxW())
0060 ,_yMax(inputParametersInstance.maxRapidity())
0061 ,_productionMode(inputParametersInstance.productionMode())
0062 ,_beamBreakupMode(inputParametersInstance.beamBreakupMode())
0063 ,_interferenceEnabled(inputParametersInstance.interferenceEnabled())
0064 ,_interferenceStrength(inputParametersInstance.interferenceStrength())
0065 ,_maxPtInterference(inputParametersInstance.maxPtInterference())
0066 ,_nmbPtBinsInterference(inputParametersInstance.nmbPtBinsInterference())
0067 ,_xsecCalcMethod(inputParametersInstance.xsecCalcMethod())
0068 ,_baseFileName(inputParametersInstance.baseFileName())
0069 {
0070   //Lets check to see if we need to recalculate the luminosity tables
0071   twoPhotonDifferentialLuminosity();
0072 }
0073
0074
0075 //______________________________________________________________________________
0076 twoPhotonLuminosity::~twoPhotonLuminosity()
0077 { }
0078
0079
0080 //______________________________________________________________________________
0081 void twoPhotonLuminosity::twoPhotonDifferentialLuminosity()
0082 {
0083   std::string wyFileName;
0084   wyFileName = _baseFileName +".txt";
0085
0086   ofstream wylumfile;
0087   wylumfile.precision(15);
0088   wylumfile.open(wyFileName.c_str());
0089   std::vector<double> w(_nWbins);
0090   std::vector<double> y(_nYbins);
0091   double xlum = 0.;
0092   double Normalize = 0.,OldNorm;
0093   double wmev = 0;
0094
0095   Normalize = 1./sqrt(1*(double)_nWbins*_nYbins); //if your grid is very fine, you'll want high accuracy-->small Normalize
0096   OldNorm   = Normalize;
0097
0098   //Writing out our input parameters+(w,y)grid+diff._lum.
0099   wylumfile << electronBeam().Z() <<endl;
0100   wylumfile << electronBeam().A() <<endl;
0101   wylumfile << targetBeam().Z() <<endl;
0102   wylumfile << targetBeam().A() <<endl;
0103   wylumfile << _gamma <<endl;
0104   wylumfile << _wMax <<endl;
0105   wylumfile << _wMin <<endl;
0106   wylumfile << _nWbins <<endl;
0107   wylumfile << _yMax <<endl;
0108   wylumfile << _nYbins <<endl;
0109   wylumfile << _productionMode <<endl;
0110   wylumfile << _beamBreakupMode <<endl;
0111   wylumfile << _interferenceEnabled <<endl;
0112   wylumfile << _interferenceStrength <<endl;
0113   wylumfile << starlightConstants::deuteronSlopePar <<endl;
0114   wylumfile << _maxPtInterference <<endl;
0115   wylumfile << _nmbPtBinsInterference <<endl;
0116   for (unsigned int i = 0; i < _nWbins; ++i) {
0117     w[i] = _wMin + (_wMax-_wMin)/_nWbins*i;
0118     wylumfile << w[i] <<endl;
0119   }
0120   for (unsigned int i = 0; i < _nYbins; ++i) {
0121     y[i] = -_yMax + 2.*_yMax*i/(_nYbins);
0122     wylumfile << y[i] <<endl;
0123   }
0124
0125   if(_xsecCalcMethod == 0) {
0126
0127     for (unsigned int i = 0; i < _nWbins; ++i) {   //For each (w,y) pair, calculate the diff. _lum
0128       for (unsigned int j = 0; j < _nYbins; ++j) {
0129         wmev = w[i]*1000.;
0130         xlum = wmev * D2LDMDY(wmev,y[j],Normalize);   //Convert photon flux dN/dW to Lorentz invariant photon number WdN/dW
0131         if (j==0) OldNorm = Normalize;       //Save value of integral for each new W(i) and Y(i)
0132         wylumfile << xlum <<endl;
0133       }
0134       Normalize = OldNorm;
0135     }
0136
0137   }
0138   else if(_xsecCalcMethod == 1) {
0139
0140     for (unsigned int i = 0; i < _nWbins; i++) {   //For each (w,y) pair, calculate the diff. _lum
0141       printf("Calculating cross section: %2.0f %% \r", float(i)/float(_nWbins)*100);
0142       fflush(stdout);
0143       for (unsigned int j = 0; j < _nYbins; j++) {
0144         xlum = w[i] * D2LDMDY(w[i],y[j]);   //Convert photon flux dN/dW to Lorentz invariant photon number WdN/dW
0145         wylumfile << xlum <<endl;
0146         // cout<<" i: "<<i<<" j: "<<j<<" W*dN/dW: "<<xlum<<endl;
0147       }
0148     }
0149   }
0150     return;
0151 }
0152
0153 //______________________________________________________________________________
0154 double twoPhotonLuminosity::D2LDMDY(double M,double Y,double &Normalize)
0155 {
0156   // double differential luminosity
0157
0158   double D2LDMDYx = 0.;
0159
0160   _W1    =  M/2.0*exp(Y);
0161   _W2    =  M/2.0*exp(-Y);
0162   int Zin1=electronBeam().Z();
0163   int Zin2=targetBeam().Z();
0164   D2LDMDYx = 2.0/M*Zin1*Zin1*Zin2*Zin2*(starlightConstants::alpha*starlightConstants::alpha)*integral(Normalize);
0165   Normalize = D2LDMDYx*M/(2.0*Zin1*Zin1*Zin2*Zin2*
0166                           starlightConstants::alpha*starlightConstants::alpha);
0167   return D2LDMDYx;
0168 }
0169
0170
0171
0172 //______________________________________________________________________________
0173 double twoPhotonLuminosity::D2LDMDY(double M, double Y) const
0174 {
0175   // double differential luminosity
0176
0177   double D2LDMDYx = 0.;
0178   double w1    =  M/2.0*exp(Y);
0179   double w2    =  M/2.0*exp(-Y);
0180
0181   double r_nuc1 = electronBeam().nuclearRadius();
0182   double r_nuc2 = targetBeam().nuclearRadius();
0183
0184   double b1min = r_nuc1;
0185   double b2min = r_nuc2;
0186
0187   double b1max = max(5.*_gamma*hbarc/w1,5*r_nuc1);
0188   double b2max = max(5.*_gamma*hbarc/w2,5*r_nuc2);
0189
0190   const int nbins_b1 = 120;
0191   const int nbins_b2 = 120;
0192
0193   double log_delta_b1 = (log(b1max)-log(b1min))/nbins_b1;
0194   double log_delta_b2 = (log(b2max)-log(b2min))/nbins_b2;
0195   double sum = 0;
0196   for(int i = 0; i < nbins_b1; ++i)
0197   {
0198       // Sum from nested integral
0199       double sum_b2 = 0;
0200       double b1_low = b1min*exp(i*log_delta_b1);
0201       double b1_high = b1min*exp((i+1)*log_delta_b1);
0202       double b1_cent = (b1_high+b1_low)/2.;
0203       for(int j = 0; j < nbins_b2; ++j)
0204       {
0205     // Sum from nested
0206     double sum_phi = 0;
0207     double b2_low = b2min*exp(j*log_delta_b2);
0208     double b2_high = b2min*exp((j+1)*log_delta_b2);
0209     double b2_cent = (b2_high+b2_low)/2.;
0210
0211     // Gaussian integration n = 10
0212     // Since cos is symmetric around 0 we only need 5 of the
0213     // points in the gaussian integration.
0214     const int ngi = 5;
0215     double weights[ngi] =
0216     {
0217       0.2955242247147529,
0218       0.2692667193099963,
0219       0.2190863625159820,
0220       0.1494513491505806,
0221       0.0666713443086881,
0222     };
0223
0224     double abscissas[ngi] =
0225     {
0226       -0.1488743389816312,
0227       -0.4333953941292472,
0228       -0.6794095682990244,
0229       -0.8650633666889845,
0230       -0.9739065285171717,
0231     };
0232
0233     for(int k = 0; k < ngi; ++k)
0234     {
0235         double b_rel = sqrt(b1_cent*b1_cent+b2_cent*b2_cent + 2.*b1_cent*b2_cent*cos(pi*(abscissas[k]+1)));
0236
0237         sum_phi += weights[k] * probabilityOfBreakup(b_rel)*2;
0238     }
0239     sum_b2 += targetBeam().photonDensity(b2_cent,w2)*pi*sum_phi*b2_cent*(b2_high-b2_low);
0240       }
0241
0242       sum += electronBeam().photonDensity(b1_cent, w1)*sum_b2*b1_cent*(b1_high-b1_low);
0243
0244   }
0245   D2LDMDYx = 2.*pi*M/2.*sum;
0246   return D2LDMDYx;
0247 }
0248
0249
0250 void * twoPhotonLuminosity::D2LDMDY_Threaded(void * a)
0251 {
0252   difflumiargs *args = (difflumiargs*)a;
0253   double M = args->m;
0254   double Y = args->y;
0255   args->res = args->self->D2LDMDY(M, Y);
0256
0257   return NULL;
0258 }
0259
0260
0261 //______________________________________________________________________________
0262 double twoPhotonLuminosity::integral(double Normalize)
0263 {
0264   int NIter = 0;
0265   int NIterMin = 0;
0266   double EPS = 0.;
0267   // double RM = 0.;
0268   double RM1 = 0.;
0269   double RM2 = 0.;
0270   double u1 = 0.;
0271   double u2 = 0.;
0272   double B1 = 0.;
0273   double B2 = 0.;
0274   double Integrala = 0.;
0275   double totsummary = 0.;
0276   double NEval      = 0.;
0277   double Lower[3];
0278   double Upper[3];
0279   double *WK = new double[500000];
0280   double Result, Summary, ResErr, NFNEVL;
0281
0282   EPS = .01*Normalize;   //This is EPS for integration, 1% of previous integral value.
0283   // Change this to the Woods-Saxon radius to be consistent with the older calculations (JN 230710)
0284   RM1  = electronBeam().nuclearRadius()/starlightConstants::hbarcmev;
0285   RM2  = targetBeam().nuclearRadius()/starlightConstants::hbarcmev;
0286   // RM  = electronBeam().woodSaxonRadius()/starlightConstants::hbarcmev;
0287
0288   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%
0289   NIterMin = 600;
0290   u1 = 9.*_gamma/_W1; //upper boundary in B1
0291   u2 = 9.*_gamma/_W2; //upper boundary in B2
0292   B1 = .4*_gamma/_W1; //intermediate boundary in B1
0293   B2 = .4*_gamma/_W2; //intermediate boundary in B2
0294   //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
0295
0296   if (u1 < RM1){
0297     Integrala = 0;
0298     totsummary = 0;
0299     NEval      = 0;
0300   }
0301   else if (B1 > RM1){
0302     if (u2 < RM2){
0303       Integrala = 0;
0304       totsummary = 0;
0305       NEval      = 0;
0306     }
0307     else if (B2 > RM2){            //integral has 4 parts
0308       Integrala = 0;
0309       totsummary = 40000;
0310       NEval      = 0;
0311       Lower[0]   = RM1;       //1
0312       Lower[1]   = RM2;       //2
0313       Lower[2]   = 0.;       //3
0314       Upper[2]   = 2.*starlightConstants::pi;    //3
0315       Upper[0]   = B1;       //1
0316       Upper[1]   = B2;       //2
0317       radmul(3,Lower,Upper,NIterMin,NIter,EPS,WK,NIter,Result,ResErr,NFNEVL,Summary);
0318       Integrala   = Integrala + Result;
0319       totsummary = totsummary + 1000*Summary;
0320       NEval      = NEval + NFNEVL;
0321       Upper[0]   = u1;       //1
0322       Upper[1]   = B2;       //2
0323       Lower[0]   = B1;       //1
0324       Lower[1]   = RM2;       //2
0325       radmul(3,Lower,Upper,NIterMin,NIter,EPS,WK,NIter,Result,ResErr,NFNEVL,Summary);
0326       Integrala   = Integrala + Result;
0327       totsummary = totsummary + 100*Summary;
0328       NEval      = NEval + NFNEVL;
0329       Upper[0]   = B1;       //1
0330       Upper[1]   = u2;       //2
0331       Lower[0]   = RM1;       //1
0332       Lower[1]   = B2;       //2
0333       radmul(3,Lower,Upper,NIterMin,NIter,EPS,WK,NIter,Result,ResErr,NFNEVL,Summary);
0334       Integrala   = Integrala + Result;
0335       totsummary = totsummary + 100*Summary;
0336       NEval      = NEval + NFNEVL;
0337       Upper[0]   = u1;       //1
0338       Upper[1]   = u2;       //2
0339       Lower[0]   = B1;       //1
0340       Lower[1]   = B2;       //2
0341       radmul(3,Lower,Upper,NIterMin,NIter,EPS,WK,NIter,Result,ResErr,NFNEVL,Summary);
0342       Integrala   = Integrala + Result;
0343       totsummary = totsummary + Summary;
0344       NEval      = NEval + NFNEVL;
0345     }
0346     else {
0347       //integral has 2 parts, b2 integral has only 1 component
0348       Integrala   = 0;
0349       totsummary = 20000;
0350       NEval      = 0;
0351       Lower[0]   = RM1;       //1
0352       Lower[1]   = RM2;       //2
0353       Lower[2]   = 0.;       //3
0354       Upper[2]   = 2.*starlightConstants::pi;    //3
0355       Upper[0]   = B1;       //1
0356       Upper[1]   = u2;       //2
0357       radmul(3,Lower,Upper,NIterMin,NIter,EPS,WK,NIter,Result,ResErr,NFNEVL,Summary);
0358       Integrala   = Integrala + Result;
0359       totsummary = totsummary + 100*Summary;
0360       NEval      = NEval + NFNEVL;
0361       Upper[0]   = u1;       //1
0362       Lower[0]   = B1;       //1
0363       radmul(3,Lower,Upper,NIterMin,NIter,EPS,WK,NIter,Result,ResErr,NFNEVL,Summary);
0364       Integrala   = Integrala + Result;
0365       totsummary = totsummary + Summary;
0366       NEval      = NEval + NFNEVL;
0367     }
0368   }
0369   else{
0370     if (u2 < RM2 ){
0371       Integrala   = 0;
0372       totsummary = 0;
0373       NEval      = 0;
0374     }
0375     else if (B2 > RM2){
0376       //integral has 2 parts, b1 integral has only 1 component
0377       Integrala   = 0;
0378       totsummary = 20000;
0379       NEval      = 0;
0380       Lower[0]   = RM1;       //1
0381       Lower[1]   = RM2;       //2
0382       Lower[2]   = 0.;       //2
0383       Upper[2]   = 2.*starlightConstants::pi;    //3
0384       Upper[0]   = u1;       //1
0385       Upper[1]   = B2;       //2
0386       radmul(3,Lower,Upper,NIterMin,NIter,EPS,WK,NIter,Result,ResErr,NFNEVL,Summary);
0387       Integrala   = Integrala + Result;
0388       totsummary = totsummary + 100*Summary;
0389       NEval      = NEval + NFNEVL;
0390       Upper[1]   = u2;       //2
0391       Lower[1]   = B2;       //2
0392       radmul(3,Lower,Upper,NIterMin,NIter,EPS,WK,NIter,Result,ResErr,NFNEVL,Summary);
0393       Integrala   = Integrala + Result;
0394       totsummary = totsummary + Summary;
0395       NEval      = NEval + NFNEVL;
0396     }
0397     else{                 //integral has 1 part
0398       Integrala   = 0;
0399       totsummary = 10000;
0400       NEval      = 0;
0401       Lower[0]   = RM1;       //1
0402       Lower[1]   = RM2;       //2
0403       Lower[2]   = 0.;       //3
0404       Upper[2]   = 2.*starlightConstants::pi;    //3
0405       Upper[0]   = u1;       //1
0406       Upper[1]   = u2;       //2
0407       radmul(3,Lower,Upper,NIterMin,NIter,EPS,WK,NIter,Result,ResErr,NFNEVL,Summary);
0408       Integrala   = Integrala + Result;
0409       totsummary = totsummary + Summary;
0410       NEval      = NEval + NFNEVL;
0411     }
0412   }
0413   Integrala = 2*starlightConstants::pi*Integrala;
0414   delete [] WK;
0415   return Integrala;
0416 }
0417
0418 //______________________________________________________________________________
0419 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)
0420 {
0421   double wn1[14] = {       -0.193872885230909911, -0.555606360818980835,
0422                -0.876695625666819078, -1.15714067977442459,  -1.39694152314179743,
0423                -1.59609815576893754,  -1.75461057765584494,  -1.87247878880251983,
0424                -1.94970278920896201,  -1.98628257887517146,  -1.98221815780114818,
0425                -1.93750952598689219,  -1.85215668343240347,  -1.72615963013768225};
0426
0427   double wn3[14] = {        0.0518213686937966768,    0.0314992633236803330,
0428                 0.0111771579535639891, -0.00914494741655235473,  -0.0294670527866686986,
0429                 -0.0497891581567850424, -0.0701112635269013768,   -0.0904333688970177241,
0430                 -0.110755474267134071,  -0.131077579637250419,    -0.151399685007366752,
0431                 -0.171721790377483099,  -0.192043895747599447,    -0.212366001117715794};
0432
0433     double wn5[14] = {        0.871183254585174982e-01,  0.435591627292587508e-01,
0434                   0.217795813646293754e-01, 0.108897906823146873e-01,  0.544489534115734364e-02,
0435                   0.272244767057867193e-02, 0.136122383528933596e-02,  0.680611917644667955e-03,
0436                   0.340305958822333977e-03, 0.170152979411166995e-03,  0.850764897055834977e-04,
0437                   0.425382448527917472e-04, 0.212691224263958736e-04,  0.106345612131979372e-04};
0438
0439     double wpn1[14] = {       -1.33196159122085045, -2.29218106995884763,
0440                   -3.11522633744855959, -3.80109739368998611, -4.34979423868312742,
0441                   -4.76131687242798352, -5.03566529492455417, -5.17283950617283939,
0442                   -5.17283950617283939, -5.03566529492455417, -4.76131687242798352,
0443                   -4.34979423868312742, -3.80109739368998611, -3.11522633744855959};
0444
0445     double wpn3[14] = {        0.0445816186556927292, -0.0240054869684499309,
0446                    -0.0925925925925925875, -0.161179698216735251,  -0.229766803840877915,
0447                    -0.298353909465020564,  -0.366941015089163228,  -0.435528120713305891,
0448                    -0.504115226337448555,  -0.572702331961591218,  -0.641289437585733882,
0449                    -0.709876543209876532,  -0.778463648834019195,  -0.847050754458161859};
0450
0451     RESULT = 0;
0452
0453     double ABSERR = 0.;
0454     double ctr[15], wth[15], wthl[15], z[15];
0455     double R1  = 1.;
0456     double HF  = R1/2.;
0457     double xl2 = 0.358568582800318073;
0458     double xl4 = 0.948683298050513796;
0459     double xl5 = 0.688247201611685289;
0460     double w2  = 980./6561.;
0461     double w4  = 200./19683.;
0462     double wp2 = 245./486.;
0463     double wp4 = 25./729.;
0464     int j1 =0;
0465     double SUM1, SUM2, SUM3, SUM4, SUM5, RGNCMP=0., RGNVAL, RGNERR, F2, F3, DIF, DIFMAX, IDVAXN =0.;
0466     IFAIL  = 3;
0467     if (N < 2 || N > 15) return 0;
0468     if (MINPTS > MAXPTS) return 0;
0469     int IFNCLS = 0;
0470     int IDVAX0 = 0;
0471     bool LDV = false;
0472     double TWONDM = pow(2.,(double)(N));
0473     int IRGNST = 2*N+3;
0474     int IRLCLS = (int)pow(2.,(N))+2*N*(N+1)+1;
0475     int ISBRGN = IRGNST;
0476     int ISBTMP, ISBTPP;
0477     int ISBRGS = IRGNST;
0478
0479     if ( MAXPTS < IRLCLS ) return 0;
0480     for ( int j = 0; j < N; j++ ){  //10
0481       ctr[j]=(B[j] + A[j])*HF;
0482       wth[j]=(B[j] - A[j])*HF;
0483     }
0484  L20:
0485     double RGNVOL = TWONDM; //20
0486     for ( int j = 0; j < N; j++ ){            //30
0487       RGNVOL = RGNVOL*wth[j];
0488       z[j] = ctr[j];
0489     }
0490
0491     SUM1 = integrand(N,z);
0492
0493     DIFMAX = 0;
0494     SUM2   = 0;
0495     SUM3   = 0;
0496     for ( int j = 0; j < N; j++ ) {   //40
0497       z[j]=ctr[j]-xl2*wth[j];
0498       F2=integrand(N,z);
0499
0500       z[j]=ctr[j]+xl2*wth[j];
0501       F2=F2+integrand(N,z);
0502       wthl[j]=xl4*wth[j];
0503
0504       z[j]=ctr[j]-wthl[j];
0505       F3=integrand(N,z);
0506
0507       z[j]=ctr[j]+wthl[j];
0508       F3=F3+integrand(N,z);
0509       SUM2=SUM2+F2;
0510       SUM3=SUM3+F3;
0511       DIF=fabs(7.*F2-F3-12.*SUM1);
0512       DIFMAX=max(DIF,DIFMAX);
0513
0514       if ( DIFMAX == DIF) IDVAXN = j+1;
0515       z[j]=ctr[j];
0516
0517     }
0518
0519     SUM4   = 0;
0520     for ( int j = 1; j < N; j++){ //70
0521
0522     j1=j-1;
0523         for ( int k = j; k < N; k++){  //60
0524       for ( int l = 0; l < 2; l++){      //50
0525         wthl[j1]=-wthl[j1];
0526         z[j1]=ctr[j1]+wthl[j1];
0527         for ( int m = 0; m < 2; m++){   //50
0528           wthl[k]=-wthl[k];
0529           z[k]=ctr[k]+wthl[k];
0530           SUM4=SUM4+integrand(N,z);
0531         }
0532       }
0533       z[k]=ctr[k];
0534     }
0535         z[j1]=ctr[j1];
0536     }
0537
0538     SUM5  = 0;
0539
0540     for ( int j = 0; j < N; j++){             //80
0541       wthl[j]=-xl5*wth[j];
0542       z[j]=ctr[j]+wthl[j];
0543     }
0544  L90:
0545     SUM5=SUM5+integrand(N,z);   //line 90
0546
0547     for (int j = 0; j < N; j++){ //100
0548       wthl[j]=-wthl[j];
0549       z[j]=ctr[j]+wthl[j];
0550       if ( wthl[j] > 0. ) goto L90;
0551     }
0552
0553     RGNCMP = RGNVOL*(wpn1[N-2]*SUM1+wp2*SUM2+wpn3[N-2]*SUM3+wp4*SUM4);
0554     RGNVAL = wn1[N-2]*SUM1+w2*SUM2+wn3[N-2]*SUM3+w4*SUM4+wn5[N-2]*SUM5;
0555
0556     RGNVAL = RGNVOL*RGNVAL;
0557     RGNERR = fabs(RGNVAL-RGNCMP);
0558     RESULT = RESULT+RGNVAL;
0559     ABSERR = ABSERR+RGNERR;
0560     IFNCLS = IFNCLS+IRLCLS;
0561
0562
0563     if (LDV){
0564     L110:
0565
0566       ISBTMP = 2*ISBRGN;
0567
0568       if ( ISBTMP > ISBRGS ) goto L160;
0569       if ( ISBTMP < ISBRGS ){
0570     ISBTPP = ISBTMP + IRGNST;
0571
0572     if ( WK[ISBTMP-1] < WK[ISBTPP-1] ) ISBTMP = ISBTPP;
0573       }
0574
0575       if ( RGNERR >= WK[ISBTMP-1] ) goto L160;
0576       for ( int k = 0; k < IRGNST; k++){
0577     WK[ISBRGN-k-1] = WK[ISBTMP-k-1];
0578                 }
0579       ISBRGN = ISBTMP;
0580       goto L110;
0581     }
0582  L140:
0583
0584     ISBTMP = (ISBRGN/(2*IRGNST))*IRGNST;
0585
0586     if ( ISBTMP >= IRGNST && RGNERR > WK[ISBTMP-1] ){
0587       for ( int k = 0; k < IRGNST; k++){
0588     WK[ISBRGN-k-1]=WK[ISBTMP-k-1];
0589       }
0590       ISBRGN = ISBTMP;
0591       goto L140;
0592     }
0593  L160:
0594
0595     WK[ISBRGN-1] = RGNERR;
0596     WK[ISBRGN-2] = RGNVAL;
0597     WK[ISBRGN-3] = IDVAXN;
0598
0599     for ( int j = 0; j < N; j++) {
0600
0601       ISBTMP = ISBRGN-2*j-4;
0602       WK[ISBTMP]=ctr[j];
0603       WK[ISBTMP-1]=wth[j];
0604     }
0605     if (LDV) {
0606       LDV = false;
0607       ctr[IDVAX0-1]=ctr[IDVAX0-1]+2*wth[IDVAX0-1];
0608       ISBRGS = ISBRGS + IRGNST;
0609       ISBRGN = ISBRGS;
0610       goto L20;
0611     }
0612
0613     RELERR=ABSERR/fabs(RESULT);
0614
0615
0616     if ( ISBRGS + IRGNST > IWK ) IFAIL = 2;
0617     if ( IFNCLS + 2*IRLCLS > MAXPTS ) IFAIL = 1;
0618     if ( RELERR < EPS && IFNCLS >= MINPTS ) IFAIL = 0;
0619
0620     if ( IFAIL == 3 ) {
0621       LDV = true;
0622       ISBRGN = IRGNST;
0623       ABSERR = ABSERR-WK[ISBRGN-1];
0624       RESULT = RESULT-WK[ISBRGN-2];
0625       IDVAX0 = (int)WK[ISBRGN-3];
0626
0627       for ( int j = 0; j < N; j++) {
0628     ISBTMP = ISBRGN-2*j-4;
0629     ctr[j] = WK[ISBTMP];
0630     wth[j] = WK[ISBTMP-1];
0631       }
0632
0633       wth[IDVAX0-1] = HF*wth[IDVAX0-1];
0634       ctr[IDVAX0-1] = ctr[IDVAX0-1]-wth[IDVAX0-1];
0635       goto L20;
0636     }
0637     NFNEVL=IFNCLS;
0638     return 1;
0639 }
0640
0641
0642 //______________________________________________________________________________
0643 double twoPhotonLuminosity::integrand(double ,  // N (unused)
0644                                       double X[15])
0645 {
0646   double  b1 = X[0];      //1
0647   double  b2 = X[1];      //2
0648   double  theta = X[2];   //3
0649   //breakup effects distances in fermis, so convert to fermis(factor of hbarcmev)
0650   double  D = sqrt(b1*b1+b2*b2-2*b1*b2*cos(theta))*starlightConstants::hbarcmev;
0651   double  integrandx = Nphoton(_W1,_gamma,b1)*Nphoton(_W2,_gamma,b2)*b1*b2*probabilityOfBreakup(D);
0652   return integrandx;
0653 }
0654
0655
0656 //______________________________________________________________________________
0657 double twoPhotonLuminosity::Nphoton(double W,double gamma,double Rho)
0658 {
0659  double Nphoton1 =0.;
0660  double WGamma = W/gamma;
0661  double WGR    = 1.0*WGamma*Rho;
0662  //factor of Z^2*alpha is omitted
0663  double Wphib = WGamma*bessel::dbesk1(WGR);
0664
0665  Nphoton1 = 1.0/(starlightConstants::pi*starlightConstants::pi)*(Wphib*Wphib);
0666  return Nphoton1;
0667 }