File indexing completed on 2024-09-27 07:03:42
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
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
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);
0096 OldNorm = Normalize;
0097
0098
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) {
0128 for (unsigned int j = 0; j < _nYbins; ++j) {
0129 wmev = w[i]*1000.;
0130 xlum = wmev * D2LDMDY(wmev,y[j],Normalize);
0131 if (j==0) OldNorm = Normalize;
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++) {
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]);
0145 wylumfile << xlum <<endl;
0146
0147 }
0148 }
0149 }
0150 return;
0151 }
0152
0153
0154 double twoPhotonLuminosity::D2LDMDY(double M,double Y,double &Normalize)
0155 {
0156
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
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
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
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
0212
0213
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
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;
0283
0284 RM1 = electronBeam().nuclearRadius()/starlightConstants::hbarcmev;
0285 RM2 = targetBeam().nuclearRadius()/starlightConstants::hbarcmev;
0286
0287
0288 NIter = 10000 + (int)1000000*(int)Normalize;
0289 NIterMin = 600;
0290 u1 = 9.*_gamma/_W1;
0291 u2 = 9.*_gamma/_W2;
0292 B1 = .4*_gamma/_W1;
0293 B2 = .4*_gamma/_W2;
0294
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){
0308 Integrala = 0;
0309 totsummary = 40000;
0310 NEval = 0;
0311 Lower[0] = RM1;
0312 Lower[1] = RM2;
0313 Lower[2] = 0.;
0314 Upper[2] = 2.*starlightConstants::pi;
0315 Upper[0] = B1;
0316 Upper[1] = B2;
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;
0322 Upper[1] = B2;
0323 Lower[0] = B1;
0324 Lower[1] = RM2;
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;
0330 Upper[1] = u2;
0331 Lower[0] = RM1;
0332 Lower[1] = B2;
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;
0338 Upper[1] = u2;
0339 Lower[0] = B1;
0340 Lower[1] = B2;
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
0348 Integrala = 0;
0349 totsummary = 20000;
0350 NEval = 0;
0351 Lower[0] = RM1;
0352 Lower[1] = RM2;
0353 Lower[2] = 0.;
0354 Upper[2] = 2.*starlightConstants::pi;
0355 Upper[0] = B1;
0356 Upper[1] = u2;
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;
0362 Lower[0] = B1;
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
0377 Integrala = 0;
0378 totsummary = 20000;
0379 NEval = 0;
0380 Lower[0] = RM1;
0381 Lower[1] = RM2;
0382 Lower[2] = 0.;
0383 Upper[2] = 2.*starlightConstants::pi;
0384 Upper[0] = u1;
0385 Upper[1] = B2;
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;
0391 Lower[1] = B2;
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{
0398 Integrala = 0;
0399 totsummary = 10000;
0400 NEval = 0;
0401 Lower[0] = RM1;
0402 Lower[1] = RM2;
0403 Lower[2] = 0.;
0404 Upper[2] = 2.*starlightConstants::pi;
0405 Upper[0] = u1;
0406 Upper[1] = u2;
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++ ){
0481 ctr[j]=(B[j] + A[j])*HF;
0482 wth[j]=(B[j] - A[j])*HF;
0483 }
0484 L20:
0485 double RGNVOL = TWONDM;
0486 for ( int j = 0; j < N; j++ ){
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++ ) {
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++){
0521
0522 j1=j-1;
0523 for ( int k = j; k < N; k++){
0524 for ( int l = 0; l < 2; l++){
0525 wthl[j1]=-wthl[j1];
0526 z[j1]=ctr[j1]+wthl[j1];
0527 for ( int m = 0; m < 2; m++){
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++){
0541 wthl[j]=-xl5*wth[j];
0542 z[j]=ctr[j]+wthl[j];
0543 }
0544 L90:
0545 SUM5=SUM5+integrand(N,z);
0546
0547 for (int j = 0; j < N; j++){
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 ,
0644 double X[15])
0645 {
0646 double b1 = X[0];
0647 double b2 = X[1];
0648 double theta = X[2];
0649
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
0663 double Wphib = WGamma*bessel::dbesk1(WGR);
0664
0665 Nphoton1 = 1.0/(starlightConstants::pi*starlightConstants::pi)*(Wphib*Wphib);
0666 return Nphoton1;
0667 }