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