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