ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Reggeons.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4Reggeons.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
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 G4ParticleDefinition * particle)
38 {
39  // Kaidalov&Piskunova Orig
40  Alpha_pomeron = 1.12; //0.9808;
41  Alphaprime_pomeron = 0.22/GeV/GeV; //0.25/GeV/GeV;
42  S0_pomeron = 1.0*GeV*GeV; //2.7*GeV*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_pomeron_NN);
54 
55  G4double Gamma_pomeron_NN = 2.14/GeV/GeV; //(2.6+3.96)
56  G4double Gamma_pomeron_N = std::sqrt(Gamma_pomeron_NN);
57  G4double Gamma_pomeron_Pr(0.), Gamma_pomeron_Tr(0.);
58 
59  G4double Rsquare_pomeron_NN = 3.30/GeV/GeV; //3.56
60  G4double Rsquare_pomeron_N = Rsquare_pomeron_NN/2.;
61  G4double Rsquare_pomeron_Pr(0.), Rsquare_pomeron_Tr(0.);
62  //-------------------------------------------------------
63 
64  // Copied from G4HadronNucleonXsc::HyperonNucleonXscNS(...)
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::HyperonNucleonXscNS(...)
76 
77  // Copied from G4HadronNucleonXsc::SCBMesonNucleonXscNS(...)
78  static const G4double llMesCof1C = 0.676568;
79  static const G4double llMesCof1B = 0.610989;
80  static const G4double llMesCof2C = 0.353135;
81  static const G4double llMesCof2B = 0.221978;
82  static const G4double llMesCofSC = 0.496568;
83  static const G4double llMesCofSB = 0.430989;
84  static const G4double llMesCofCB = 0.287557;
85  static const G4double llMesCofEtaP = 0.88;
86  static const G4double llMesCofEta = 0.76;
87  // End of Copied from G4HadronNucleonXsc::SCBMesonNucleonXscNS(...)
88 
89  if ( absPDGcode > 1000 ) { // Projectile is baryon or anti_baryon --------
90 
91  Cpr_pomeron = C_pomeron_N; // Shower enhancement coefficient for projectile
92  Ctr_pomeron = C_pomeron_N; // Shower enhancement coefficient for target
93  C_pomeron = Cpr_pomeron*Ctr_pomeron;
94 
95  Gamma_pomeron_Pr = Gamma_pomeron_N; // vertex constant for projectile
96  Gamma_pomeron_Tr = Gamma_pomeron_N; // vertex constant for target
97  Gamma_pomeron = Gamma_pomeron_Pr * Gamma_pomeron_Tr;
98 
99  Rsquare_pomeron_Pr = Rsquare_pomeron_N; // R^2 of pomeron-projectile interaction
100  Rsquare_pomeron_Tr = Rsquare_pomeron_N; // R^2 of pomeron-target interaction
101  Rsquare_pomeron = Rsquare_pomeron_Pr + Rsquare_pomeron_Tr;
102 
103  Freggeon_Alpha = 0.7; // Intersept of f-trajectory
104  Freggeon_Alphaprime = 0.8/GeV/GeV; // Slope of f-trajectory
105  Freggeon_Gamma = sqr(2.871)/GeV/GeV; // Vertex constant of f-meson - nucleon interactions
106  Freggeon_Rsquare = 2*0.916/GeV/GeV; // R^2 of f-meson - nucleon interactions
107  Freggeon_C = 1.0; // Shower enhancement coefficient
108  FParity = +1; // Parity of the trajectory
109 
110  Wreggeon_Alpha = 0.4; // Intersept of omega-trajectory (w)
111  Wreggeon_Alphaprime = 0.9/GeV/GeV; // Slope of w-trajectory
112  Wreggeon_Gamma = sqr(2.241)/GeV/GeV; // Vertex constant of w-meson - nucleon interactions
113  Wreggeon_Rsquare = 2*0.945/GeV/GeV *0.5; // R^2 of w-meson - nucleon interactions
114  Wreggeon_C = 1.0; // Shower enhancement coefficient
115  if (PDGcode > 0) WParity = -1; // Parity +1 for Pbar P, and -1 for PP interactions
116  if (PDGcode < 0) WParity = +1;
117 
118  // Copied from G4HadronNucleonXsc::HyperonNucleonXscNS(...)
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 == 4432 || PDGcode ==-4432 ) {
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::HyperonNucleonXscNS(...)
159 
160  Gamma_pomeron_Pr *= coeff;
161 
162  } else if ( absPDGcode == 211 || PDGcode == 111 || absPDGcode >= 400 ) { // Projectile is a nonstrange meson
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 -> 0.89 // Uzhi
169  Gamma_pomeron_Tr = Gamma_pomeron_N;
170  Gamma_pomeron = Gamma_pomeron_Pr * Gamma_pomeron_Tr;
171 
172  Rsquare_pomeron_Pr = 0.5/GeV/GeV;
173  Rsquare_pomeron_Tr = Rsquare_pomeron_N;
174  Rsquare_pomeron = Rsquare_pomeron_Pr + Rsquare_pomeron_Tr;
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; // 1.12 -> 0.56 Uzhi
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::SCBMesonNucleonXscNS(...)
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::SCBMesonNucleonXscNS(...)
222 
223  Gamma_pomeron_Pr *= coeff;
224 
225  } else if ( absPDGcode == 321 || absPDGcode == 311 ||
226  PDGcode == 130 || PDGcode == 310 ) { // Projectile is a Kaon
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 * Gamma_pomeron_Tr;
235 
236  Rsquare_pomeron_Pr = 0.31/GeV/GeV;
237  Rsquare_pomeron_Tr = Rsquare_pomeron_N;
238  Rsquare_pomeron = Rsquare_pomeron_Pr + Rsquare_pomeron_Tr;
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 ) { // Projectile is Gamma
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; // 1.415
264  Gamma_pomeron_Tr = Gamma_pomeron_N;
265  Gamma_pomeron = Gamma_pomeron_Pr * Gamma_pomeron_Tr;
266 
267  Rsquare_pomeron_Pr = 0.51/GeV/GeV;
268  Rsquare_pomeron_Tr = Rsquare_pomeron_N;
269  Rsquare_pomeron = Rsquare_pomeron_Pr + Rsquare_pomeron_Tr;
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, Nucleon assumed
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 * Gamma_pomeron_Tr;
294 
295  Rsquare_pomeron_Pr = Rsquare_pomeron_N;
296  Rsquare_pomeron_Tr = Rsquare_pomeron_N;
297  Rsquare_pomeron = Rsquare_pomeron_Pr + Rsquare_pomeron_Tr;
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.; XDDiff=0.;
318  Xinel =0.; Xnd=0.; XndP=0.; XndR=0.;
319  /*
320  G4cout<<G4endl<<"Reggeon's parameters for Particle "<<particle->GetParticleName()<<" "<<PDGcode<<G4endl<<G4endl;
321  G4cout<<"Alpha_pomeron "<<Alpha_pomeron;
322  G4cout<<" Alphaprime_pomeron "<<Alphaprime_pomeron*GeV*GeV;
323  G4cout<<" S0_pomeron "<<S0_pomeron/GeV/GeV<<G4endl;
324  G4cout<<"Gamma_pomeron "<<Gamma_pomeron*GeV*GeV;
325  G4cout<<" Rsquare_pomeron "<<Rsquare_pomeron*GeV*GeV;
326  G4cout<<" C_pomeron "<<C_pomeron<<G4endl<<G4endl;
327  G4int Uzhi; G4cin>>Uzhi;
328  */
329 }
330 
331 
332 G4double G4Reggeons::Get_Cprojectile() {return Cpr_pomeron;}
333 
334 
335 G4double G4Reggeons::Get_Ctarget() {return Ctr_pomeron;}
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.; XDDiff=0.; G4double XDiff=0.;
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.); // chiPin Pomeron inelastic phase is a data member
358  chiPin=0.;
359  for(G4int i=0; i<10000;i++)
360  {
361  B += dB;
362 
363  chiP = Chi_pomeron(1.,B); chiR = Chi_reggeon(1.,B);
364  chiPin = Chi_pomeron(2.,B); chiRin = Chi_reggeon(2.,B);
365 
366  AmplitudeP = (1.0/C_pomeron)*(1.0 - G4Exp(-chiP))*G4Exp(-chiR);
367  AmplitudeR = (1.0 - G4Exp(-chiR));
368 
369  Xtotal += 2 * (AmplitudeP + AmplitudeR) * B * dB;
370  XtotalP += 2 * (AmplitudeP + 0. ) * B * dB;
371  XtotalR += 2 * (0. + AmplitudeR) * B * dB;
372 
373  Xelastic += sqr(AmplitudeP + AmplitudeR) * B * dB;
374  Xpr_Diff += (Cpr_pomeron - 1.0) * sqr(AmplitudeP) * B * dB;
375  Xtr_Diff += (Ctr_pomeron - 1.0) * sqr(AmplitudeP) * B * dB;
376  XDiff += (Cpr_pomeron - 1.0) * (Ctr_pomeron - 1.0) * sqr(AmplitudeP) * B * dB;
377 
378  // ----------------------------------
379  AmplitudeP = (1.0/C_pomeron)*(1.0 - G4Exp(-chiPin))*G4Exp(-chiRin);
380  AmplitudeR = (1.0 - G4Exp(-chiRin));
381 
382  Xnd += (AmplitudeP + AmplitudeR) * B * dB;
383  XndP += (AmplitudeP + 0. ) * B * dB;
384  XndR += (0. + AmplitudeR) * B * dB;
385  }
386 
387  Xtotal *=twopi; XtotalP *=twopi; XtotalR *=twopi;
388  Xelastic *=twopi; Xpr_Diff *=twopi; Xtr_Diff *=twopi; XDiff *=twopi;
389  Xinel = Xtotal - Xelastic;
390  (void)Xinel; // To avoid compiler warning "variable not used"
391 
392  Xnd *=twopi; XndP *=twopi; XndR *=twopi;
393  XDDiff = XDiff-Xpr_Diff-Xtr_Diff;
394 
395  /*
396  G4cout<<"Total totalP totalR "<<Xtotal/millibarn <<" "<<XtotalP/millibarn <<" "<<XtotalR/millibarn<<" mb"<<G4endl;
397  G4cout<<"Elastic "<<Xelastic/millibarn <<G4endl;
398  G4cout<<"PrDiff TrDiff W_Diff "<<Xpr_Diff/millibarn<<" "<<Xtr_Diff/millibarn<<" "<<XDiff/millibarn<<G4endl;
399  G4cout<<"Inelastic "<<Xinel/millibarn <<G4endl;
400  G4cout<<"NonDiff Pom & Reg "<<Xnd/millibarn <<" "<<XndP/millibarn <<" "<<XndR/millibarn <<G4endl;
401  */
402 }
403 
404 
405 G4double G4Reggeons::Chi_pomeron(G4double Mult, G4double B)
406 {
407  G4double R2 = Rsquare_pomeron + Alphaprime_pomeron * G4Log(Sint/S0_pomeron);
408  G4double Eikonal = Mult * C_pomeron * Gamma_pomeron/R2 *
409  G4Pow::GetInstance()->powA(Sint/S0_pomeron, Alpha_pomeron -1.) *
410  G4Exp(-sqr(B)/4.0/R2/hbarc_squared);
411  return Eikonal;
412 }
413 
414 
415 G4double G4Reggeons::Chi_reggeon(G4double Mult, G4double B)
416 {
417  G4double R2F = Freggeon_Rsquare + Freggeon_Alphaprime * G4Log(Sint/S0_pomeron);
418  G4double R2W = Wreggeon_Rsquare + Wreggeon_Alphaprime * G4Log(Sint/S0_pomeron);
419 
420  G4double Eikonal = Mult * FParity * Freggeon_C * Freggeon_Gamma/R2F *
421  G4Pow::GetInstance()->powA(Sint/S0_pomeron, Freggeon_Alpha -1.) *
422  G4Exp(-sqr(B)/4.0/R2F/hbarc_squared);
423 
424  Eikonal+= Mult * WParity * Wreggeon_C * Wreggeon_Gamma/R2W *
425  G4Pow::GetInstance()->powA(Sint/S0_pomeron, Wreggeon_Alpha -1.) *
426  G4Exp(-sqr(B)/4.0/R2W/hbarc_squared);
427  return Eikonal;
428 }
429 
430 
431 G4double G4Reggeons::GetTotalX() { return Xtotal; }
432 G4double G4Reggeons::GetTotalXp() { return XtotalP; }
433 G4double G4Reggeons::GetTotalXr() { return XtotalR; }
434 
435 G4double G4Reggeons::GetElasticX(){ return Xelastic; }
436 G4double G4Reggeons::GetPrDiffX() { return Xpr_Diff; }
437 G4double G4Reggeons::GetTrDiffX() { return Xtr_Diff; }
438 G4double G4Reggeons::GetDDiffX() { return XDDiff; }
439 
440 G4double G4Reggeons::GetInelX() { return Xinel; }
441 
442 G4double G4Reggeons::GetND_X() { return Xnd; }
443 G4double G4Reggeons::GetNDp_X() { return XndP; }
444 G4double G4Reggeons::GetNDr_X() { return XndR; }
445 
446 
447 //----------------------------------------------------------------------------------------------
448 void G4Reggeons::GetProbabilities(G4double B, G4int Mode,
449  G4double & Pint,
450  G4double & Pprd, G4double & Ptrd, G4double & Pdd,
451  G4double & Pnd, G4double & Pnvr)
452 {
453  // Puprose of the method is a calculation of inelastic interaction probability (Pint),
454  // probability of projectile diffraction (Pprd),
455  // probability of target diffraction (Ptrd),
456  // probability of double diffraction (Pdd ),
457  // probability of non-diffractive inter. (Pnd ),
458  // probability of quark-exc. inter. (Pnvr),
459  // number of cutted pomerons (NcutPomerons).
460  // The input parameters are B - impact parameter, and Mode = All/WITHOUT_R/NON_DIFF
461  //
462  if ( B > 2.* fermi ) { Pint=0.; Pprd=0.; Ptrd=0.; Pdd=0.; Pnd=0.; Pnvr=0.; return;}
463  // At large B for hN collisions it is better to return zero inter. probability
464 
465  G4double chiP = Chi_pomeron(1.,B); G4double chiR = Chi_reggeon(1.,B);
466  chiPin = Chi_pomeron(2.,B); G4double chiRin = Chi_reggeon(2.,B);
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 - G4Exp(-chiP))*Exp_ChiR;
472  G4double AmplitudeR = (1.0 - Exp_ChiR);
473 
474  G4double AmplitudeP2, Apr_Diff, Atr_Diff, ADiff;
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_pomeron - 1.0) * AmplitudeP2;
481 
482  // ----------------------------------
483  Exp_ChiR = G4Exp(-chiRin);
484  AmplitudeP = (1.0/C_pomeron)*(1.0 - G4Exp(-chiPin))*Exp_ChiR;
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; // Probability of projectile diffraction
498  Ptrd = Atr_Diff/Pint; // Probability of target diffraction
499  Pdd = ADiff /Pint; // Probability of double diffraction
500  Pnd = AndP /Pint; // Probability of non-diffractive inelastic interaction
501  Pnvr = AndR /Pint; // Probability of non-vacuum reggeon (nvr) inelastic interaction
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-complete Poisson distribution
526 {
527  if ( chiPin < 0.001 ) return 0; // At small average multiplicity of cutted pomerons
528  // it is better to return 0 to avoid problems with
529  // calculation exactness.
530  G4double ksi = G4UniformRand() * (1.0-G4Exp(-chiPin)) * G4Exp(chiPin);
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