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