]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STARLIGHT/starlight/src/beambeamsystem.cpp
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update to trunk of hepforge
[u/mrichter/AliRoot.git] / STARLIGHT / starlight / src / beambeamsystem.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:
45d54d9a 23// $Rev:: 176 $: revision of last commit
24// $Author:: jseger $: author of last commit
25// $Date:: 2014-06-20 22:15:20 +0200 #$: date of last commit
da32329d
AM		 26//
27// Description:
28//
29//
30//
31///////////////////////////////////////////////////////////////////////////
32
33
34#include <iostream>
35#include <fstream>
36#include <cmath>
37
38#include "inputParameters.h"
39#include "reportingUtils.h"
40#include "starlightconstants.h"
41#include "bessel.h"
42#include "beambeamsystem.h"
43
44
45using namespace std;
46using namespace starlightConstants;
47
48
49//______________________________________________________________________________
45d54d9a 50beamBeamSystem::beamBeamSystem(const beam& beam1,
51 const beam& beam2)
da32329d
AM		 52 : _beamLorentzGamma(inputParametersInstance.beamLorentzGamma()),
53 _beamBreakupMode (inputParametersInstance.beamBreakupMode()),
45d54d9a 54 _beam1 (beam1),
55 _beam2 (beam2),
da32329d
AM		 56 _breakupProbabilities(0),
57 _breakupImpactParameterStep(1.007),
58 _breakupCutOff(10e-6)
59{
60 init();
61}
62
63
64
65
66//______________________________________________________________________________
67beamBeamSystem::beamBeamSystem()
68 : _beamLorentzGamma(inputParametersInstance.beamLorentzGamma()),
69 _beamBreakupMode (inputParametersInstance.beamBreakupMode()),
70 _beam1 (inputParametersInstance.beam1Z(),
71 inputParametersInstance.beam1A(),
da32329d
AM		 72 inputParametersInstance.coherentProduction()),
73 _beam2 (inputParametersInstance.beam2Z(),
74 inputParametersInstance.beam2A(),
da32329d
AM		 75 inputParametersInstance.coherentProduction()),
76 _breakupProbabilities(0),
77 _breakupImpactParameterStep(1.007),
78 _breakupCutOff(10e-10)
79{
80 init();
81}
82
83
84
85//______________________________________________________________________________
86beamBeamSystem::~beamBeamSystem()
87{ }
88
89void beamBeamSystem::init()
90{
91 // Calculate beam gamma in CMS frame
92 double rap1 = acosh(inputParametersInstance.beam1LorentzGamma());
93 double rap2 = -acosh(inputParametersInstance.beam2LorentzGamma());
94
95 _cmsBoost = (rap1+rap2)/2.;
96
97 _beamLorentzGamma = cosh((rap1-rap2)/2);
98 _beam1.setBeamLorentzGamma(_beamLorentzGamma);
99 _beam2.setBeamLorentzGamma(_beamLorentzGamma);
100
101 generateBreakupProbabilities();
102}
103//______________________________________________________________________________
104double
105beamBeamSystem::probabilityOfBreakup(const double D) const
106{
107
108 double bMin = (_beam1.nuclearRadius()+_beam2.nuclearRadius())/2.;
109 double pOfB = 0.; // PofB = 1 means that there will be a UPC event and PofB = 0 means no UPC
110
111 // Do pp here
112 if ((_beam1.Z() == 1) && (_beam1.A() == 1) && (_beam2.Z() == 1) && (_beam2.A() == 1)) {
113 double ppslope=19.8;
114 double GammaProfile = exp(-D * D / (2. * hbarc * hbarc * ppslope));
115 pOfB = (1. - GammaProfile) * (1. - GammaProfile);
116 // if (D < 2. * _beam1.nuclearRadius())
117 // //Should be the total of RNuc1+Rnuc2,used only beam #1
118 // PofB = 0.0;
119 // else
120 // PofB = 1.0;
121 return pOfB;
122 }
123 else if ( ( (_beam1.A() == 1) && (_beam2.A() != 1) ) || ((_beam1.A() != 1) && (_beam2.A() == 1)) ) {
124 // This is pA
125 if( _beam1.A() == 1 ){
45d54d9a 126 bMin = _beam2.nuclearRadius() + 0.7;
127 pOfB = exp(-7.0*_beam2.rho0()*_beam2.thickness(D));
da32329d
AM		 128 }else if( _beam2.A() == 1 ){
129 bMin = _beam1.nuclearRadius() + 0.7;
45d54d9a 130 pOfB = exp(-7.0*_beam1.rho0()*_beam1.thickness(D));
da32329d
AM		 131 }else{
132 cout<<"Some logical problem here!"<<endl;
133 }
45d54d9a 134 // if( D > bMin )pOfB=1.0;
da32329d 135 return pOfB;
45d54d9a 136
da32329d
AM		 137 }
138
139 //use the lookup table and return...
140 pOfB = 1.;
141 if (D > 0.0) {
142 //Now we must determine which step number in d corresponds to this D,
143 // and use appropiate Ptot(D_i)
144 //int i = (int)(log(D / Bmin) / log(1.01));
145 int i = (int)(log(D / bMin) / log(_breakupImpactParameterStep));
146 if (i <= 0)
147 pOfB = _breakupProbabilities[0];
148 else{
149 if (i >= int(_breakupProbabilities.size()-1))
150 pOfB = _breakupProbabilities[_breakupProbabilities.size()-1];
151 else {
152 // const double DLow = Bmin * pow((1.01), i);
153 const double DLow = bMin * pow((_breakupImpactParameterStep), i);
154 // const double DeltaD = 0.01 * DLow;
155 const double DeltaD = (_breakupImpactParameterStep-1) * DLow;
156 const double DeltaP = _breakupProbabilities[i + 1] - _breakupProbabilities[i];
157 pOfB = _breakupProbabilities[i] + DeltaP * (D - DLow) / DeltaD;
158 }
159 }
160 }
161
162 return pOfB;
163}
164
165void
166beamBeamSystem::generateBreakupProbabilities()
167{
168 // Step = 1.007;//.01; //We will multiplicateively increase Biter by 1%
169
170
171 double bMin = (_beam1.nuclearRadius()+_beam2.nuclearRadius())/2.;
172
173
174 if ((_beam1.Z() != 1) && (_beam1.A() != 1) && (_beam2.Z() != 1) && _beam2.A() != 1) {
175
176 if (_beamBreakupMode == 1)
177 printInfo << "Hard Sphere Break criteria. b > " << 2. * _beam1.nuclearRadius() << endl;
178 if (_beamBreakupMode == 2)
179 printInfo << "Requiring XnXn [Coulomb] breakup. " << endl;
180 if (_beamBreakupMode == 3)
181 printInfo << "Requiring 1n1n [Coulomb only] breakup. " << endl;
182 if (_beamBreakupMode == 4)
183 printInfo << "Requiring both nuclei to remain intact. " << endl;
184 if (_beamBreakupMode == 5)
185 printInfo << "Requiring no hadronic interactions. " << endl;
186 if (_beamBreakupMode == 6)
187 printInfo << "Requiring breakup of one or both nuclei. " << endl;
188 if (_beamBreakupMode == 7)
189 printInfo << "Requiring breakup of one nucleus (Xn,0n). " << endl;
190
191 //pp may cause segmentation fault in here and it does not use this...
192 double pOfB = 0;
193 double b = bMin;
194 double totRad = _beam1.nuclearRadius()+_beam2.nuclearRadius();
195
196 while(1)
197 {
198
199 if(_beamBreakupMode != 5)
200 {
201 if(b > (totRad*1.5))
202 {
203 if(pOfB<_breakupCutOff)
204 {
205// std::cout << "Break off b: " << b << std::endl;
206// std::cout << "Number of PofB bins: " << _breakupProbabilities.size() << std::endl;
207 break;
208 }
209 }
210 }
211 else
212 {
213 if((1-pOfB)<_breakupCutOff)
214 {
215 // std::cout << "Break off b: " << b << std::endl;
216// std::cout << "Number of PofB bins: " << _breakupProbabilities.size() << std::endl;
217 break;
218 }
219 }
220// std::cout << 1-pOfBreakup << std::endl;
221// _pHadronBreakup = 0;
222// _pPhotonBreakup = 0;
223
224// double pHadronBreakup = probabilityOfHadronBreakup(b);
225 probabilityOfHadronBreakup(b);
226 //moved gammatarg into photonbreakup
227// double pPhotonBreakup = probabilityOfPhotonBreakup(b, _beamBreakupMode);
228 probabilityOfPhotonBreakup(b, _beamBreakupMode);
229
230 //What was probability of photonbreakup depending upon mode selection,
231 // is now done in the photonbreakupfunction
232 if (_beamBreakupMode == 1) {
233 if (b >_beam1.nuclearRadius()+_beam2.nuclearRadius()) // symmetry
234 _pHadronBreakup = 0;
235 else
236 _pHadronBreakup = 999.;
237 }
238
239 b *= _breakupImpactParameterStep;
240 pOfB = exp(-1 * _pHadronBreakup) * _pPhotonBreakup;
241 _breakupProbabilities.push_back(pOfB);
242 }
243 }
244 else if (((_beam1.Z() == 1) && (_beam1.A() == 1)) || ((_beam2.Z() == 1) && (_beam2.A() == 1))) {
245
246 double pOfB = 0;
247 double b = bMin;
248 double totRad = _beam1.nuclearRadius()+_beam2.nuclearRadius();
249
250 while(1)
251 {
252 if(_beamBreakupMode != 5)
253 {
254 if(b > (totRad*1.5))
255 {
256 if(pOfB<_breakupCutOff)
257 {
258// std::cout << "Break off b: " << b << std::endl;
259 break;
260 }
261 }
262 }
263 else
264 {
265 if((1-pOfB)<_breakupCutOff)
266 {
267// std::cout << "Break off b: " << b << std::endl;
268 break;
269 }
270 }
271 _beam1.Z() > 1 ? pOfB = exp(-7.35*_beam1.thickness(b)) :
272 pOfB = exp(-7.35*_beam2.thickness(b));
273 _breakupProbabilities.push_back(pOfB);
274 b *= _breakupImpactParameterStep;
275 }
276 }
277
278
279}
280
281//______________________________________________________________________________
282double
283beamBeamSystem::probabilityOfHadronBreakup(const double impactparameter)
284{
285 // double pbreakup =0.;
286 //probability of hadron breakup,
287 //this is what is returned when the function is called
288 double gamma = _beamLorentzGamma;
289 //input for gamma_em
290 //this will need to be StarlightInputParameters::gamma or whatever
291 double b = impactparameter;
292 int ap = _beam1.A();
293 //Notice this is taking from nucleus 1.Still assuming symmetric system?
294
295 static int IFIRSTH = 0;
296 static double DELL=0., DELR=0., SIGNN=0., R1=0., A1=0., R2=0., RHO1=0.;
297 static double RHO2=0., NZ1=0., NZ2=0., NR1=0., NR2=0.,RR1=0., NY=0., NX=0.;
298 static double AN1=0., AN2=0.;
299 double delo=0.,RSQ=0.,Z1=0.,Y=0.,X=0.,XB=0.,RPU=0.,IRUP=0.,RTU=0.;
300 double IRUT=0.,T1=0.,T2=0.;
301 static double DEN1[20002], DEN2[20002];
302 if (IFIRSTH != 0) goto L100;
303 //Initialize
304 //Integration delta x, delta z
305 IFIRSTH = 1;
306 DELL = .05;
307 DELR = .01;
308
309 //use two sigma_NN's. 52mb at rhic 100gev/beam, 88mb at LHC 2.9tev/beam, gamma is in cm system
310 SIGNN = 5.2;
311 if ( gamma > 500. ) SIGNN = 8.8;
312 //use parameter from Constants
313 R1 = ( _beam1.nuclearRadius()); //remember _beam2? better way to do this generically
314 A1 = 0.535; //This is woodsaxonskindepth?
315 //write(6,12)r1,a1,signn Here is where we could probably set this up asymmetrically R2=_beam2.nuclearRadius() and RHO2=ap2=_beam2.A()
316 R2 = R1;
317 RHO1 = ap;
318 RHO2 = RHO1;
319 NZ1 = ((R1+5.)/DELR);
320 NR1 = NZ1;
321 NZ2 = ((R2+5.)/DELR);
322 NR2 = NZ2;
323 RR1 = -DELR;
324 NY = ((R1+5.)/DELL);
325 NX = 2*NY;
326 for ( int IR1 = 1; IR1 <= NR1; IR1++) {
327 DEN1[IR1] = 0.;
328 RR1 = RR1+DELR;
329 Z1 = -DELR/2;
330
331 for ( int IZ1 = 1; IZ1 <= NZ1; IZ1++) {
332 Z1 = Z1+DELR;
333 RSQ = RR1*RR1+Z1*Z1;
334 DEN1[IR1] = DEN1[IR1]+1./(1.+exp((sqrt(RSQ)-R1)/A1));
335 }//for(IZ)
336
337 DEN1[IR1] = DEN1[IR1]*2.*DELR;
338 DEN2[IR1] = DEN1[IR1];
339 }//for(IR)
340
341 AN1 = 0.;
342 RR1 = 0.;
343
344 for ( int IR1 =1; IR1 <= NR1; IR1++) {
345 RR1 = RR1+DELR;
346 AN1 = AN1+RR1*DEN1[IR1]*DELR*2.*3.141592654;
347 }//for(IR)
348
349 AN2 = AN1; //This will also probably need to be changed?
350
351 delo = .05;
352 //.1 to turn mb into fm^2
353 //Calculate breakup probability here
354 L100:
355 _pHadronBreakup = 0.;
356 if ( b > 25. ) return _pHadronBreakup;
357 Y = -.5*DELL;
358 for ( int IY = 1; IY <= NY; IY++) {
359 Y = Y+DELL;
360 X = -DELL*float(NY+1);
361
362 for ( int IX = 1; IX <=NX; IX++) {
363 X = X+DELL;
364 XB = b-X;
365 RPU = sqrt(X*X+Y*Y);
366 IRUP = (RPU/DELR)+1;
367 RTU = sqrt(XB*XB+Y*Y);
368 IRUT = (RTU/DELR)+1;
369 T1 = DEN2[(int)IRUT]*RHO2/AN2;
370 T2 = DEN1[(int)IRUP]*RHO1/AN1;
371 //Eq.6 BCW, Baltz, Chasman, White, Nucl. Inst. & Methods A 417, 1 (1998)
372 _pHadronBreakup=_pHadronBreakup+2.*T1*(1.-exp(-SIGNN*T2))*DELL*DELL;
373 }//for(IX)
374 }//for(IY)
375
376 return _pHadronBreakup;
377}
378
379
380//______________________________________________________________________________
381double
382beamBeamSystem::probabilityOfPhotonBreakup(const double impactparameter, const int mode)
383{
384 static double ee[10001], eee[162], se[10001];
385
386 //double gamma_em=108.4; //This will be an input value.
387 _pPhotonBreakup =0.; //Might default the probability with a different value?
388 double b = impactparameter;
389 int zp = _beam1.Z(); //What about _beam2? Generic approach?
390 int ap = _beam1.A();
391
392 //Was initialized at the start of the function originally, been moved inward.
393 double pxn=0.;
394 double p1n=0.;
395
396 //Used to be done prior to entering the function. Done properly for assymetric?
397 double gammatarg = 2.*_beamLorentzGamma*_beamLorentzGamma-1.;
398 double omaxx =0.;
399 //This was done prior entering the function as well
400 if (_beamLorentzGamma > 500.){
401 omaxx=1.E10;
402 }
403 else{
404 omaxx=1.E7;
405 }
406
407
408 double e1[23]= {0.,103.,106.,112.,119.,127.,132.,145.,171.,199.,230.,235.,
409 254.,280.,300.,320.,330.,333.,373.,390.,420.,426.,440.};
410 double s1[23]= {0.,12.0,11.5,12.0,12.0,12.0,15.0,17.0,28.0,33.0,
411 52.0,60.0,70.0,76.0,85.0,86.0,89.0,89.0,75.0,76.0,69.0,59.0,61.0};
412 double e2[12]={0.,2000.,3270.,4100.,4810.,6210.,6600.,
413 7790.,8400.,9510.,13600.,16400.};
414 double s2[12]={0.,.1266,.1080,.0805,.1017,.0942,.0844,.0841,.0755,.0827,
415 .0626,.0740};
416 double e3[29]={0.,26.,28.,30.,32.,34.,36.,38.,40.,44.,46.,48.,50.,52.,55.,
417 57.,62.,64.,66.,69.,72.,74.,76.,79.,82.,86.,92.,98.,103.};
418 double s3[29]={0.,30.,21.5,22.5,18.5,17.5,15.,14.5,19.,17.5,16.,14.,
419 20.,16.5,17.5,17.,15.5,18.,15.5,15.5,15.,13.5,18.,14.5,15.5,12.5,13.,
420 13.,12.};
421 static double sa[161]={0.,0.,.004,.008,.013,.017,.021,.025,.029,.034,.038,.042,.046,
422 .051,.055,.059,.063,.067,.072,.076,.08,.085,.09,.095,.1,.108,.116,
423 .124,.132,.14,.152,.164,.176,.188,.2,.22,.24,.26,.28,.3,.32,.34,
424 .36,.38,.4,.417,.433,.450,.467,.483,.5,.51,.516,.52,.523,.5245,
425 .525,.5242,
426 .5214,.518,.512,.505,.495,.482,.469,.456,.442,.428,.414,.4,.386,
427 .370,.355,.34,.325,.310,.295,.280,.265,.25,.236,.222,.208,.194,
428 .180,.166,
429 .152,.138,.124,.11,.101,.095,.09,.085,.08,.076,.072,.069,.066,
430 .063,.06,.0575,.055,.0525,.05,.04875,.0475,.04625,.045,.04375,
431 .0425,.04125,.04,.03875,.0375,.03625,.035,.03375,.0325,.03125,.03,
432 .02925,.0285,.02775,.027,.02625,.0255,.02475,.024,.02325,.0225,
433 .02175,.021,.02025,.0195,.01875,.018,.01725,.0165,.01575,.015,
434 .01425,.0135,.01275,.012,.01125,.0105,.00975,.009,.00825,.0075,
435 .00675,.006,.00525,.0045,.00375,.003,.00225,.0015,.00075,0.};
436
437
438
439 double sen[161]={0.,0.,.012,.025,.038,.028,.028,.038,.035,.029,.039,.035,
440 .038,.032,.038,.041,.041,.049,.055,.061,.072,.076,.070,.067,
441 .080,.103,.125,.138,.118,.103,.129,.155,.170,.180,.190,.200,
442 .215,.250,.302,.310,.301,.315,.330,.355,.380,.400,.410,.420,
443 .438,.456,.474,.492,.510,.533,.556,.578,.6,.62,.63,.638,
444 .640,.640,.637,.631,.625,.618,.610,.600,.580,.555,.530,.505,
445 .480,.455,.435,.410,.385,.360,.340,.320,.300,.285,.270,.255,
446 .240,.225,.210,.180,.165,.150,.140,.132,.124,.116,.108,.100,
447 .092,.084,.077,.071,.066,.060,.055,.051,.048,.046,.044,.042,
448 .040,.038,.036,.034,.032,.030,.028,.027,.026,.025,.025,.025,
449 .024,.024,.024,.024,.024,.023,.023,.023,.023,.023,.022,.022,
450 .022,.022,.022,.021,.021,.021,.020,.020,
451 .020,.019,.018,.017,.016,.015,.014,.013,.012,.011,.010,.009,
452 .008,.007,.006,.005,.004,.003,.002,.001,0.};
453
454 // gammay,p gamma,n of Armstrong begin at 265 incr 25
455
456
457 double sigt[160]={0.,.4245,.4870,.5269,.4778,.4066,.3341,.2444,.2245,.2005,
458 .1783,.1769,.1869,.1940,.2117,.2226,.2327,.2395,.2646,.2790,.2756,
459 .2607,.2447,.2211,.2063,.2137,.2088,.2017,.2050,.2015,.2121,.2175,
460 .2152,.1917,.1911,.1747,.1650,.1587,.1622,.1496,.1486,.1438,.1556,
461 .1468,.1536,.1544,.1536,.1468,.1535,.1442,.1515,.1559,.1541,.1461,
462 .1388,.1565,.1502,.1503,.1454,.1389,.1445,.1425,.1415,.1424,.1432,
463 .1486,.1539,.1354,.1480,.1443,.1435,.1491,.1435,.1380,.1317,.1445,
464 .1375,.1449,.1359,.1383,.1390,.1361,.1286,.1359,.1395,.1327,.1387,
465 .1431,.1403,.1404,.1389,.1410,.1304,.1363,.1241,.1284,.1299,.1325,
466 .1343,.1387,.1328,.1444,.1334,.1362,.1302,.1338,.1339,.1304,.1314,
467 .1287,.1404,.1383,.1292,.1436,.1280,.1326,.1321,.1268,.1278,.1243,
468 .1239,.1271,.1213,.1338,.1287,.1343,.1231,.1317,.1214,.1370,.1232,
469 .1301,.1348,.1294,.1278,.1227,.1218,.1198,.1193,.1342,.1323,.1248,
470 .1220,.1139,.1271,.1224,.1347,.1249,.1163,.1362,.1236,.1462,.1356,
471 .1198,.1419,.1324,.1288,.1336,.1335,.1266};
472
473
474 double sigtn[160]={0.,.3125,.3930,.4401,.4582,.3774,.3329,.2996,.2715,.2165,
475 .2297,.1861,.1551,.2020,.2073,.2064,.2193,.2275,.2384,.2150,.2494,
476 .2133,.2023,.1969,.1797,.1693,.1642,.1463,.1280,.1555,.1489,.1435,
477 .1398,.1573,.1479,.1493,.1417,.1403,.1258,.1354,.1394,.1420,.1364,
478 .1325,.1455,.1326,.1397,.1286,.1260,.1314,.1378,.1353,.1264,.1471,
479 .1650,.1311,.1261,.1348,.1277,.1518,.1297,.1452,.1453,.1598,.1323,
480 .1234,.1212,.1333,.1434,.1380,.1330,.12,.12,.12,.12,.12,.12,.12,.12,
481 .12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,
482 .12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,
483 .12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,
484 .12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,
485 .12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12};
486 //89*.12};
487
488
489
490 static int IFIRSTP=0;
491
492
493 // Initialization needed?
494 //double hbar=197.3;
495 //double pi=3.141592654;
496
497 // added
498 double si1=0, g1 =0, o1=0;
499 int ne = 0, ij =0;
500 double delo=0, omax =0, gk1m=0;
501 static double scon=0., zcon=0.,o0=0.;
502
503
504 double x=0,y=0,eps=0,eta=0,em=0,exx=0,s=0,ictr=0,pom=0,vec=0,gk1=0;
505
506 // maximum energy for GDR dissocation (in target frame, in MeV)
507
508 double omax1n=24.01;
509
510 if (IFIRSTP != 0) goto L100;
511
512 IFIRSTP=1;
513
514
515 //This whole thing is dependenant on gold or lead....Might need to expand
516 if (zp == 79)
517 {
518
519
520 ap=197;
521 si1=540.;
522 g1=4.75;
523
524 // peak and minimum energies for GDR excitation (in MeV)
525 o1=13.70;
526 o0=8.1;
527 }
528 else
529 {
530 zp=82; //assumed to be lead
531 ap=208;
532 si1=640.;
533 g1=4.05;
534 o1=13.42;
535 o0=7.4;
536 for(int j=1;j<=160;j++)
537 {
538
539 sa[j]=sen[j];
540 }
541 }
542 //Part II of initialization
543 delo = .05;
544 //.1 to turn mb into fm^2
545 scon = .1*g1*g1*si1;
546 zcon = zp/(gammatarg*( pi)*(
547 hbarcmev))*zp/(gammatarg*( pi)*
548 ( hbarcmev))/137.04;
549
550 //single neutron from GDR, Veyssiere et al. Nucl. Phys. A159, 561 (1970)
551 for ( int i = 1; i <= 160; i++) {
552 eee[i] = o0+.1*(i-1);
553 sa[i] = 100.*sa[i];
554 }
555 //See Baltz, Rhoades-Brown, and Weneser, Phys. Rev. E 54, 4233 (1996)
556 //for details of the following photo cross-sections
557 eee[161]=24.1;
558 ne=int((25.-o0)/delo)+1;
559 //GDR any number of neutrons, Veyssiere et al., Nucl. Phys. A159, 561 (1970)
560 for ( int i = 1; i <= ne; i++ ) {
561 ee[i] = o0+(i-1)*delo;
562 //cout<<" ee 1 "<<ee[i]<<" "<<i<<endl;
563
564 se[i] = scon*ee[i]*ee[i]/(((o1*o1-ee[i]*ee[i])*(o1*o1-ee[i]*ee[i]))
565 +ee[i]*ee[i]*g1*g1);
566 }
567 ij = ne; //Risky?
568 //25-103 MeV, Lepretre, et al., Nucl. Phys. A367, 237 (1981)
569 for ( int j = 1; j <= 27; j++ ) {
570 ij = ij+1;
571 ee[ij] = e3[j];
572 //cout<<" ee 2 "<<ee[ij]<<" "<<ij<<endl;
573
574 se[ij] = .1*ap*s3[j]/208.;
575 }
576 //103-440 MeV, Carlos, et al., Nucl. Phys. A431, 573 (1984)
577 for ( int j = 1; j <= 22; j++ ) {
578 ij = ij+1;
579 ee[ij] = e1[j];
580 //cout<<" ee 3 "<<ee[ij]<<" "<<ij<<endl;
581 se[ij] = .1*ap*s1[j]/208.;
582 }
583 //440 MeV-2 GeV Armstrong et al.
584 for ( int j = 9; j <= 70; j++) {
585 ij = ij+1;
586 ee[ij] = ee[ij-1]+25.;
587 //cout<<" ee 4 "<<ee[ij]<<" "<<ij<<endl;
588 se[ij] = .1*(zp*sigt[j]+(ap-zp)*sigtn[j]);
589 }
590 //2-16.4 GeV Michalowski; Caldwell
591 for ( int j = 1; j <= 11; j++) {
592 ij = ij+1;
593 ee[ij] = e2[j];
594 //cout<<" ee 5 "<<ee[ij]<<" "<<ij<<endl;
595 se[ij] = .1*ap*s2[j];
596 }
597 //Regge paramteres
598 x = .0677;
599 y = .129;
600 eps = .0808;
601 eta = .4525;
602 em = .94;
603 exx = pow(10,.05);
604
605 //Regge model for high energy
606 s = .002*em*ee[ij];
607 //make sure we reach LHC energies
608 ictr = 100;
609 if ( gammatarg > (2.*150.*150.)) ictr = 150;
610 for ( int j = 1; j <= ictr; j++ ) {
611 ij = ij+1;
612 s = s*exx;
613 ee[ij] = 1000.*.5*(s-em*em)/em;
614 //cout<<" ee 6 "<<ee[ij]<<" "<<ij<<endl;
615 pom = x*pow(s,eps);
616 vec = y*pow(s,(-eta));
617 se[ij] = .1*.65*ap*(pom+vec);
618 }
619 ee[ij+1] = 99999999999.;
620 //done with initaliation
621 //write(6,99)o0;
622 //clear counters for 1N, XN
623 L100:
624
625 p1n = 0.;
626 pxn = 0.;
627 //start XN calculation
628 //what's the b-dependent highest energy of interest?
629
630 omax = min(omaxx,4.*gammatarg*( hbarcmev)/b);
631 if ( omax < o0 ) return _pPhotonBreakup;
632 gk1m = bessel::dbesk1(ee[1]*b/(( hbarcmev)*gammatarg));
633 int k = 2;
634 L212:
635 if (ee[k] < omax ) {
636 gk1 = bessel::dbesk1(ee[k]*b/(( hbarcmev)*gammatarg));
637 //Eq. 3 of BCW--NIM in Physics Research A 417 (1998) pp1-8:
638 pxn=pxn+zcon*(ee[k]-ee[k-1])*.5*(se[k-1]*ee[k-1]*gk1m*gk1m+se[k]*ee[k]*gk1*gk1);
639 k = k + 1;
640 gk1m = gk1;
641 goto L212;
642 }
643 //one neutron dissociation
644 omax = min(omax1n,4.*gammatarg*( hbarcmev)/b);
645 gk1m = bessel::dbesk1(eee[1]*b/(( hbarcmev)*gammatarg));
646 k = 2;
647 L102:
648 if (eee[k] < omax ) {
649 gk1 = bessel::dbesk1(eee[k]*b/(( hbarcmev)*gammatarg));
650 //Like Eq3 but with only the one neutron out GDR photo cross section input
651 p1n = p1n+zcon*(eee[k]-eee[k-1])*.5*(sa[k-1]*eee[k-1]*gk1m*gk1m+sa[k]*eee[k]*gk1*gk1);
652 k = k+1;
653 gk1m = gk1;
654 goto L102;
655 }
656
657
658 //This used to be done externally, now it is done internally.
659 if (( mode) == 1) _pPhotonBreakup = 1.;
660 if (( mode) == 2) _pPhotonBreakup = (1-exp(-1*pxn))*(1-exp(-1*pxn));
661 if (( mode) == 3) _pPhotonBreakup = (p1n*exp(-1*pxn))*(p1n*exp(-1*pxn));
662 if (( mode) == 4) _pPhotonBreakup = exp(-2*pxn);
663 if (( mode) == 5) _pPhotonBreakup = 1.;
664 if (( mode) == 6) _pPhotonBreakup = (1. - exp(-2.*pxn));
665 if (( mode) == 7) _pPhotonBreakup = 2.*exp(-pxn)*(1.-exp(-pxn));
666
667 //cout<<pxn<<" "<<zcon<<" "<<ee[k]<<" "<<se[k-1]<<" "<<gk1m<<" "<<gk1<<" "<<k<<" "<<ee[k+1]<< " "<<b<< endl;
668
669 return _pPhotonBreakup;
670}