| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 //
28
29 #include "G4Reggeons.hh"
30 #include "G4PhysicalConstants.hh"
31 #include "G4SystemOfUnits.hh"
32 #include "G4Pow.hh"
33 #include "G4Exp.hh"
34 #include "G4Log.hh"
35
36
37 G4Reggeons::G4Reggeons(const G4ParticleDefinit
38 {
39 // Kaidalov&Piskunova Orig
40 Alpha_pomeron = 1.12; //0.9808;
41 Alphaprime_pomeron = 0.22/GeV/GeV; //0.25/Ge
42 S0_pomeron = 1.0*GeV*GeV; //2.7*GeV
43
44 Alpha_pomeronHard = 1.47;
45 Gamma_pomeronHard = 0.0/GeV/GeV;
46
47 G4int PDGcode = particle->GetPDGEncoding();
48 G4int absPDGcode = std::abs(PDGcode);
49
50 //------------------------------------------
51 // Kaidalov&Piskunova Orig
52 G4double C_pomeron_NN = 1.5;
53 G4double C_pomeron_N = std::sqrt(C_pome
54
55 G4double Gamma_pomeron_NN = 2.14/GeV/GeV;
56 G4double Gamma_pomeron_N = std::sqrt(Gamma_
57 G4double Gamma_pomeron_Pr(0.), Gamma_pomeron
58
59 G4double Rsquare_pomeron_NN = 3.30/GeV/GeV;
60 G4double Rsquare_pomeron_N = Rsquare_pomero
61 G4double Rsquare_pomeron_Pr(0.), Rsquare_pom
62 //------------------------------------------
63
64 // Copied from G4HadronNucleonXsc::HyperonNu
65 G4double coeff = 1.0;
66 static const G4double lBarCof1S = 0.88;
67 static const G4double lBarCof2S = 0.76;
68 static const G4double lBarCof3S = 0.64;
69 static const G4double lBarCof1C = 0.784378
70 static const G4double lBarCofSC = 0.664378
71 static const G4double lBarCof2SC = 0.544378
72 static const G4double lBarCof1B = 0.740659
73 static const G4double lBarCofSB = 0.620659
74 static const G4double lBarCof2SB = 0.500659
75 // End of Copied from G4HadronNucleonXsc::Hy
76
77 // Copied from G4HadronNucleonXsc::SCBMesonN
78 static const G4double llMesCof1C = 0.67656
79 static const G4double llMesCof1B = 0.61098
80 static const G4double llMesCof2C = 0.35313
81 static const G4double llMesCof2B = 0.22197
82 static const G4double llMesCofSC = 0.49656
83 static const G4double llMesCofSB = 0.43098
84 static const G4double llMesCofCB = 0.28755
85 static const G4double llMesCofEtaP = 0.88;
86 static const G4double llMesCofEta = 0.76;
87 // End of Copied from G4HadronNucleonXsc::SC
88
89 if ( absPDGcode > 1000 ) { // Projectile i
90
91 Cpr_pomeron = C_pomeron_N;
92 Ctr_pomeron = C_pomeron_N;
93 C_pomeron = Cpr_pomeron*Ctr_pomeron;
94
95 Gamma_pomeron_Pr = Gamma_pomeron_N;
96 Gamma_pomeron_Tr = Gamma_pomeron_N;
97 Gamma_pomeron = Gamma_pomeron_Pr * Gamm
98
99 Rsquare_pomeron_Pr = Rsquare_pomeron_N;
100 Rsquare_pomeron_Tr = Rsquare_pomeron_N;
101 Rsquare_pomeron = Rsquare_pomeron_Pr +
102
103 Freggeon_Alpha = 0.7;
104 Freggeon_Alphaprime = 0.8/GeV/GeV;
105 Freggeon_Gamma = sqr(2.871)/GeV/GeV
106 Freggeon_Rsquare = 2*0.916/GeV/GeV;
107 Freggeon_C = 1.0;
108 FParity = +1;
109
110 Wreggeon_Alpha = 0.4;
111 Wreggeon_Alphaprime = 0.9/GeV/GeV;
112 Wreggeon_Gamma = sqr(2.241)/GeV/GeV;
113 Wreggeon_Rsquare = 2*0.945/GeV/GeV *0.
114 Wreggeon_C = 1.0;
115 if (PDGcode > 0) WParity = -1;
116 if (PDGcode < 0) WParity = +1;
117
118 // Copied from G4HadronNucleonXsc::Hyperon
119 if ( PDGcode == 3122 || PDGcode == 3222 ||
120 PDGcode == 3112 || PDGcode == 3212 ||
121 PDGcode ==-3122 || PDGcode ==-3222 ||
122 PDGcode ==-3112 || PDGcode ==-3212
123 coeff = lBarCof1S;
124 }
125 if ( PDGcode == 3312 || PDGcode == 3322 ||
126 PDGcode ==-3312 || PDGcode ==-3322
127 coeff = lBarCof2S;
128 }
129 if ( PDGcode == 3334 || PDGcode ==-3334 )
130 coeff = lBarCof3S;
131 }
132 if ( PDGcode == 4122 || PDGcode ==-4122 ||
133 PDGcode == 4222 || PDGcode ==-4222 || /
134 PDGcode == 4212 || PDGcode ==-4212 || /
135 PDGcode == 4112 || PDGcode ==-4112 ) { /
136 coeff = lBarCof1C;
137 }
138 if ( PDGcode == 4332 || PDGcode ==-4332 )
139 coeff = lBarCof2SC;
140 }
141 if ( PDGcode == 4232 || PDGcode == 4132 ||
142 PDGcode ==-4232 || PDGcode ==-4132
143 coeff = lBarCofSC;
144 }
145 if ( PDGcode == 5122 || PDGcode ==-5122 ||
146 PDGcode == 5222 || PDGcode ==-5222 || /
147 PDGcode == 5112 || PDGcode ==-5112 || /
148 PDGcode == 5212 || PDGcode ==-5212 ) { /
149 coeff = lBarCof1B;
150 }
151 if ( PDGcode == 5332 || PDGcode ==-5332 )
152 coeff = lBarCof2SB;
153 }
154 if ( PDGcode == 5132 || PDGcode == 5232 ||
155 PDGcode ==-5132 || PDGcode ==-5232
156 coeff = lBarCofSB;
157 }
158 // End of Copied from G4HadronNucleonXsc::
159
160 Gamma_pomeron_Pr *= coeff;
161
162 } else if ( absPDGcode == 211 || PDGcode ==
163
164 Cpr_pomeron = 1.352;
165 Ctr_pomeron = C_pomeron_N;
166 C_pomeron = Cpr_pomeron*Ctr_pomeron;
167 // KP
168 Gamma_pomeron_Pr = 0.89/GeV; // 0.85 ->
169 Gamma_pomeron_Tr = Gamma_pomeron_N;
170 Gamma_pomeron = Gamma_pomeron_Pr * Gamm
171
172 Rsquare_pomeron_Pr = 0.5/GeV/GeV;
173 Rsquare_pomeron_Tr = Rsquare_pomeron_N;
174 Rsquare_pomeron = Rsquare_pomeron_Pr +
175
176 Freggeon_Alpha = 0.7;
177 Freggeon_Alphaprime = 0.8/GeV/GeV;
178 Freggeon_Gamma = 3.524/GeV/GeV;
179 Freggeon_Rsquare = 1.0/GeV/GeV;
180 Freggeon_C = 1.0;
181 FParity = +1;
182
183 Wreggeon_Alpha = 0.5;
184 Wreggeon_Gamma = 0.56/GeV/GeV; //
185 Wreggeon_Rsquare = 9.19/GeV/GeV;
186 Wreggeon_Alphaprime = 0.9/GeV/GeV;
187 Wreggeon_C = 1.0;
188 if (PDGcode > 0) WParity = -1;
189 if (PDGcode < 0) WParity = +1;
190
191 // Copied from G4HadronNucleonXsc::SCBMeso
192 if ( PDGcode == 511 || PDGcode ==-511 ||
193 PDGcode == 521 || PDGcode ==-521 )
194 coeff = llMesCof1B;
195 }
196 if ( PDGcode == 421 || PDGcode ==-421 ||
197 PDGcode == 411 || PDGcode ==-411 )
198 coeff = llMesCof1C;
199 }
200 if ( PDGcode == 531 || PDGcode ==-531 ) {
201 coeff = llMesCofSB;
202 }
203 if ( PDGcode == 541 || PDGcode ==-541 ) {
204 coeff = llMesCofCB;
205 }
206 if ( PDGcode == 431 || PDGcode ==-431 ) {
207 coeff = llMesCofSC;
208 }
209 if ( PDGcode == 441 || PDGcode == 443 ) {
210 coeff = llMesCof2C;
211 }
212 if ( PDGcode == 553 ) {
213 coeff = llMesCof2B;
214 }
215 if ( PDGcode == 221 ) {
216 coeff = llMesCofEta;
217 }
218 if ( PDGcode == 331 ) {
219 coeff = llMesCofEtaP;
220 }
221 // End of Copied from G4HadronNucleonXsc::
222
223 Gamma_pomeron_Pr *= coeff;
224
225 } else if ( absPDGcode == 321 || absPDGcode
226 PDGcode == 130 || PDGcode
227
228 Cpr_pomeron = 1.522;
229 Ctr_pomeron = C_pomeron_N;
230 C_pomeron = Cpr_pomeron*Ctr_pomeron;
231
232 Gamma_pomeron_Pr = 0.90/GeV;
233 Gamma_pomeron_Tr = Gamma_pomeron_N;
234 Gamma_pomeron = Gamma_pomeron_Pr * Gamm
235
236 Rsquare_pomeron_Pr = 0.31/GeV/GeV;
237 Rsquare_pomeron_Tr = Rsquare_pomeron_N;
238 Rsquare_pomeron = Rsquare_pomeron_Pr +
239
240 Freggeon_Alpha = 0.7;
241 Freggeon_Alphaprime = 0.8/GeV/GeV;
242 Freggeon_Gamma = 1.32/GeV/GeV;
243 Freggeon_Rsquare = 0.5/GeV/GeV;
244
245 Freggeon_C = 1.0;
246 FParity = +1;
247
248 Wreggeon_Alpha = 0.4;
249 Wreggeon_Alphaprime = 0.9 /GeV/GeV;
250 Wreggeon_Gamma = 1.68/GeV/GeV;
251 Wreggeon_Rsquare = 9.19/GeV/GeV;
252 Wreggeon_C = 1.0;
253
254 if (PDGcode > 0) WParity = -1;
255 if (PDGcode < 0) WParity = +1;
256
257 } else if ( absPDGcode == 22 ) { // Project
258
259 Cpr_pomeron = 1.437;
260 Ctr_pomeron = C_pomeron_N;
261 C_pomeron = Cpr_pomeron*Ctr_pomeron;
262
263 Gamma_pomeron_Pr = 0.0035/GeV; //
264 Gamma_pomeron_Tr = Gamma_pomeron_N;
265 Gamma_pomeron = Gamma_pomeron_Pr * Gamm
266
267 Rsquare_pomeron_Pr = 0.51/GeV/GeV;
268 Rsquare_pomeron_Tr = Rsquare_pomeron_N;
269 Rsquare_pomeron = Rsquare_pomeron_Pr +
270
271 Freggeon_Alpha = 0.7;
272 Freggeon_Alphaprime = 0.8/GeV/GeV;
273 Freggeon_Gamma = 0.011/GeV/GeV;
274 Freggeon_Rsquare = 0.5/GeV/GeV;
275 Freggeon_C = 1.0;
276 FParity = +1;
277
278 Wreggeon_Alpha = 0.;
279 Wreggeon_Alphaprime = 0.9/GeV/GeV;
280 Wreggeon_Gamma = 0.01/GeV/GeV;
281 Wreggeon_Rsquare = 1.0/GeV/GeV;
282 Wreggeon_C = 1.0;
283 WParity = +1;
284
285 } else { // Projectile is undefined, Nucleo
286
287 Cpr_pomeron = C_pomeron_N;
288 Ctr_pomeron = C_pomeron_N;
289 C_pomeron = Cpr_pomeron*Ctr_pomeron;
290
291 Gamma_pomeron_Pr = Gamma_pomeron_N;
292 Gamma_pomeron_Tr = Gamma_pomeron_N;
293 Gamma_pomeron = Gamma_pomeron_Pr * Gamm
294
295 Rsquare_pomeron_Pr = Rsquare_pomeron_N;
296 Rsquare_pomeron_Tr = Rsquare_pomeron_N;
297 Rsquare_pomeron = Rsquare_pomeron_Pr +
298
299 Freggeon_Alpha = 0.723;
300 Freggeon_Gamma = 8.801/GeV/GeV;
301 Freggeon_Rsquare = 0.396/GeV/GeV;
302 Freggeon_Alphaprime = 1.324/GeV/GeV;
303 Freggeon_C = 1.0;
304 FParity = +1;
305
306 Wreggeon_Alpha = 0.353;
307 Wreggeon_Gamma = 8.516/GeV/GeV;
308 Wreggeon_Rsquare = 24.40/GeV/GeV;
309 Wreggeon_Alphaprime = 1.5/GeV/GeV;
310 Wreggeon_C = 1.0;
311 WParity = -1;
312
313 }
314
315 chiPin=0.;
316 Xtotal =0.; XtotalP=0.; XtotalR=0.;
317 Xelastic=0.; Xpr_Diff=0.; Xtr_Diff=0.; XDDif
318 Xinel =0.; Xnd=0.; XndP=0.; XndR=0.;
319 /*
320 G4cout<<G4endl<<"Reggeon's parameters for Pa
321 G4cout<<"Alpha_pomeron "<<Alpha_pomeron;
322 G4cout<<" Alphaprime_pomeron "<<Alphaprime_p
323 G4cout<<" S0_pomeron "<<S0_pomeron/GeV/GeV<<
324 G4cout<<"Gamma_pomeron "<<Gamma_pomeron*GeV*
325 G4cout<<" Rsquare_pomeron "<<Rsquare_pomeron
326 G4cout<<" C_pomeron "<<C_pomeron<<G4endl
327 G4int Uzhi; G4cin>>Uzhi;
328 */
329 }
330
331
332 G4double G4Reggeons::Get_Cprojectile() {return
333
334
335 G4double G4Reggeons::Get_Ctarget() {return
336
337
338 G4Reggeons::~G4Reggeons() {}
339
340
341 void G4Reggeons::SetS(G4double S) {Sint = S;}
342
343
344 void G4Reggeons::CalculateXs()
345 {
346 Xtotal =0.; XtotalP=0.; XtotalR=0.;
347 Xelastic=0.; Xpr_Diff=0.; Xtr_Diff=0.; XDDif
348 Xinel =0.; Xnd=0.; XndP=0.; XndR=0.;
349
350 G4double AmplitudeP(0.), AmplitudeR(0.);
351
352 G4double B_max = 10.*fermi;
353 G4double dB = B_max/10000.;
354
355 G4double B =-dB/2.;
356
357 G4double chiP(0.), chiR(0.), chiRin(0.); /
358 chiPin=0.;
359 for(G4int i=0; i<10000;i++)
360 {
361 B += dB;
362
363 chiP = Chi_pomeron(1.,B); chiR = Chi_r
364 chiPin = Chi_pomeron(2.,B); chiRin = Chi_r
365
366 AmplitudeP = (1.0/C_pomeron)*(1.0 - G4Exp(
367 AmplitudeR = (1.0 - G4Exp(
368
369 Xtotal += 2 * (AmplitudeP + AmplitudeR)
370 XtotalP += 2 * (AmplitudeP + 0. )
371 XtotalR += 2 * (0. + AmplitudeR)
372
373 Xelastic += sqr(Ampl
374 Xpr_Diff += (Cpr_pomeron - 1.0) * sqr(Amp
375 Xtr_Diff += (Ctr_pomeron - 1.0) * sqr(Amp
376 XDiff += (Cpr_pomeron - 1.0) * (Ctr_po
377
378 // ----------------------------------
379 AmplitudeP = (1.0/C_pomeron)*(1.0 - G4Exp(
380 AmplitudeR = (1.0 - G4Exp(
381
382 Xnd += (AmplitudeP + AmplitudeR) * B
383 XndP += (AmplitudeP + 0. ) * B
384 XndR += (0. + AmplitudeR) * B
385 }
386
387 Xtotal *=twopi; XtotalP *=twopi; XtotalR *=t
388 Xelastic *=twopi; Xpr_Diff *=twopi; Xtr_Diff
389 Xinel = Xtotal - Xelastic;
390 (void)Xinel; // To avoid compiler warning "
391
392 Xnd *=twopi; XndP *=twopi; XndR *=twopi;
393 XDDiff = XDiff-Xpr_Diff-Xtr_Diff;
394
395 /*
396 G4cout<<"Total totalP totalR "<<Xtotal/mill
397 G4cout<<"Elastic "<<Xelastic/mi
398 G4cout<<"PrDiff TrDiff W_Diff "<<Xpr_Diff/mi
399 G4cout<<"Inelastic "<<Xinel/milli
400 G4cout<<"NonDiff Pom & Reg "<<Xnd/milliba
401 */
402 }
403
404
405 G4double G4Reggeons::Chi_pomeron(G4double Mult
406 {
407 G4double R2 = Rsquare_pomeron + Alphaprime_p
408 G4double Eikonal = Mult * C_pomeron * Gamma_
409 G4Pow::GetInstance()->powA(Sint/S0_po
410 G4Exp(-sqr(B)/4.0/R2/hbar
411 return Eikonal;
412 }
413
414
415 G4double G4Reggeons::Chi_reggeon(G4double Mult
416 {
417 G4double R2F = Freggeon_Rsquare + Freggeon_A
418 G4double R2W = Wreggeon_Rsquare + Wreggeon_A
419
420 G4double Eikonal = Mult * FParity * Freggeon
421 G4Pow::GetInstance()->powA(Sint/S0_po
422 G4Exp(-sqr(B)/4.0/R2F/hba
423
424 Eikonal+= Mult * WParity * Wreggeon_C * Wreg
425 G4Pow::GetInstance()->powA(Sint/S0_pomer
426 G4Exp(-sqr(B)/4.0/R2W/hbarc_square
427 return Eikonal;
428 }
429
430
431 G4double G4Reggeons::GetTotalX() { return Xto
432 G4double G4Reggeons::GetTotalXp() { return Xto
433 G4double G4Reggeons::GetTotalXr() { return Xto
434
435 G4double G4Reggeons::GetElasticX(){ return Xel
436 G4double G4Reggeons::GetPrDiffX() { return Xpr
437 G4double G4Reggeons::GetTrDiffX() { return Xtr
438 G4double G4Reggeons::GetDDiffX() { return XDD
439
440 G4double G4Reggeons::GetInelX() { return Xin
441
442 G4double G4Reggeons::GetND_X() { return Xnd
443 G4double G4Reggeons::GetNDp_X() { return Xnd
444 G4double G4Reggeons::GetNDr_X() { return Xnd
445
446
447 //--------------------------------------------
448 void G4Reggeons::GetProbabilities(G4double B,
449 G4double & Pint,
450 G4double & P
451 G4double & P
452 {
453 // Puprose of the method is a calculation of
454 //
455 //
456 //
457 //
458 //
459 //
460 // The input parameters are B - impact param
461 //
462 if ( B > 2.* fermi ) { Pint=0.; Pprd=0.; Ptr
463 // At large B for hN collisions it is better
464
465 G4double chiP = Chi_pomeron(1.,B); G4doubl
466 chiPin = Chi_pomeron(2.,B); G4doubl
467 //chiPin is data member of the class
468
469 G4double Exp_ChiR = G4Exp(-chiR);
470
471 G4double AmplitudeP = (1.0/C_pomeron)*(1.0 -
472 G4double AmplitudeR = (1.0 -
473
474 G4double AmplitudeP2, Apr_Diff, Atr_Diff, AD
475
476 // Aelastic = sqr(AmplitudeP + AmplitudeR);
477 AmplitudeP2 = sqr(AmplitudeP);
478 Apr_Diff = (Cpr_pomeron - 1.0) * AmplitudeP2
479 Atr_Diff = (Ctr_pomeron - 1.0) * AmplitudeP2
480 ADiff = (Cpr_pomeron - 1.0) * (Ctr_pomero
481
482 // ----------------------------------
483 Exp_ChiR = G4Exp(-chiRin);
484 AmplitudeP = (1.0/C_pomeron)*(1.0 - G4Exp(-c
485 AmplitudeR = (1.0 - Exp_ChiR
486
487 G4double And, AndP, AndR;
488
489 And = (AmplitudeP + AmplitudeR);
490 AndP = (AmplitudeP + 0. );
491 AndR = (0. + AmplitudeR);
492
493 // ----------------------------------
494 if ( Mode == ALL)
495 {
496 Pint = Apr_Diff + Atr_Diff + ADiff + And;
497 Pprd = Apr_Diff/Pint;
498 Ptrd = Atr_Diff/Pint;
499 Pdd = ADiff /Pint;
500 Pnd = AndP /Pint;
501 Pnvr = AndR /Pint;
502 }
503 else if ( Mode == WITHOUT_R)
504 {
505 Pint = Apr_Diff + Atr_Diff + ADiff + AndP;
506 Pprd = Apr_Diff/Pint;
507 Ptrd = Atr_Diff/Pint;
508 Pdd = ADiff /Pint;
509 Pnd = AndP /Pint;
510 Pnvr = 0.;
511 }
512 else
513 { // Mode == NON_DIFF (of projectile)
514 Pint = Atr_Diff + AndP;
515 Pprd = 0.;
516 Ptrd = Atr_Diff/Pint;
517 Pdd = 0.;
518 Pnd = AndP /Pint;
519 Pnvr = 0.;
520 }
521 return;
522 }
523
524
525 G4int G4Reggeons::ncPomerons() // Non-
526 {
527 if ( chiPin < 0.001 ) return 0; // At small
528 // it is better to return 0 to avoi
529 // calculation exactness.
530 G4double ksi = G4UniformRand() * (1.0-G4Exp
531 G4double Term = chiPin;
532 G4double Sum = Term;
533 G4int nCuts = 1;
534
535 while ( Sum < ksi) {
536 nCuts++;
537 Term *= chiPin/(G4double) nCuts;
538 Sum += Term;
539 }
540
541 return nCuts;
542 }
543
544