]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STARLIGHT/starlight/src/twophotonluminosity.cpp
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
STARLIGHT code and interface
[u/mrichter/AliRoot.git] / STARLIGHT / starlight / src / twophotonluminosity.cpp
CommitLineData
da32329d
AM		 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:: 164 $: revision of last commit
24// $Author:: odjuvsla $: author of last commit
25// $Date:: 2013-10-06 16:18:08 +0200 #$: 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
47using namespace std;
48using namespace starlightConstants;
49
50
51
52//______________________________________________________________________________
53twoPhotonLuminosity::twoPhotonLuminosity(beam beam_1,beam beam_2):
54beamBeamSystem(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//______________________________________________________________________________
69twoPhotonLuminosity::~twoPhotonLuminosity()
70{ }
71
72
73//______________________________________________________________________________
74void 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//______________________________________________________________________________
186double 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//______________________________________________________________________________
206double 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
286void * 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//______________________________________________________________________________
308double 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//______________________________________________________________________________
461double 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//______________________________________________________________________________
685double 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//______________________________________________________________________________
701double 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}