54{
55
56
57
58
59
60
61
62
63
64
65 if( vecLen == 0 )return true;
66
67
68
69 if( currentParticle.
GetMass() == 0.0 || targetParticle.
GetMass() == 0.0 )
return true;
70
74 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
75 targetMass*targetMass +
76 2.0*targetMass*etOriginal );
78 G4double availableEnergy = centerofmassEnergy-(targetMass+currentMass);
79 if( availableEnergy <= 1.0 )return true;
80
97
100
101
102
103 const G4double avrs[] = {3.,4.,5.,6.,7.,8.,9.,10.,20.,30.,40.,50.};
104
108
111 ed <<
" While count exceeded " <<
G4endl;
112 while ((i<12) && (centerofmassEnergy>avrs[i]) ) {
113 i++;
114 loop++;
115 if (loop > 1000) {
117 break;
118 }
119 }
120
121
122 if( i == 12 )
123 ibin = 11;
124 else
125 ibin = i;
126
127
128
129
130 if (vecLen == 1) {
131 i3 = 0;
132 i4 = 1;
133 } else {
136 i4 = i3 =
G4int(vecLen*rnd);
137
138 while(i3 == i4) {
141 i4 =
G4int(vecLen*rnd);
142 }
143 }
144
145
146
147 const G4double avkkb[] = { 0.0015, 0.005, 0.012, 0.0285, 0.0525, 0.075,
148 0.0975, 0.123, 0.28, 0.398, 0.495, 0.573 };
149 const G4double avky[] = { 0.005, 0.03, 0.064, 0.095, 0.115, 0.13,
150 0.145, 0.155, 0.20, 0.205, 0.210, 0.212 };
151 const G4double avnnb[] = { 0.00001, 0.0001, 0.0006, 0.0025, 0.01, 0.02,
152 0.04, 0.05, 0.12, 0.15, 0.18, 0.20 };
153
154 avk = (
G4Log(avkkb[ibin])-
G4Log(avkkb[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
155 /(avrs[ibin]-avrs[ibin-1]) +
G4Log(avkkb[ibin-1]);
157
158 avy = (
G4Log(avky[ibin])-
G4Log(avky[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
159 /(avrs[ibin]-avrs[ibin-1]) +
G4Log(avky[ibin-1]);
161
162 avn = (
G4Log(avnnb[ibin])-
G4Log(avnnb[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
163 /(avrs[ibin]-avrs[ibin-1]) +
G4Log(avnnb[ibin-1]);
165
166 if( avk+avy+avn <= 0.0 )return true;
167
168 if( currentMass < protonMass )avy /= 2.0;
169 if( targetMass < protonMass )avy = 0.0;
170 avy += avk+avn;
171 avk += avn;
173 if( ran < avn )
174 {
175 if( availableEnergy < 2.0 )return true;
176 if( vecLen == 1 )
177 {
180 {
181 vec[0]->SetDefinition( aNeutron );
184 vec[0]->SetMayBeKilled(false);
186 }
187 else
188 {
189 vec[0]->SetDefinition( aProton );
192 vec[0]->SetMayBeKilled(false);
194 }
196
197 }
198 else
199 {
201 {
202 vec[i3]->SetDefinition( aNeutron );
203 vec[i4]->SetDefinition( anAntiNeutron );
204 vec[i3]->SetMayBeKilled(false);
205 vec[i4]->SetMayBeKilled(false);
206 }
207 else
208 {
209 vec[i3]->SetDefinition( aProton );
210 vec[i4]->SetDefinition( anAntiProton );
211 vec[i3]->SetMayBeKilled(false);
212 vec[i4]->SetMayBeKilled(false);
213 }
214 }
215 }
216 else if( ran < avk )
217 {
218 if( availableEnergy < 1.0 )return true;
219
220 const G4double kkb[] = { 0.2500, 0.3750, 0.5000, 0.5625, 0.6250,
221 0.6875, 0.7500, 0.8750, 1.000 };
222 const G4int ipakkb1[] = { 10, 10, 10, 11, 11, 12, 12, 11, 12 };
223 const G4int ipakkb2[] = { 13, 11, 12, 11, 12, 11, 12, 13, 13 };
225 i = 0;
226
227 loop = 0;
229 eda <<
" While count exceeded " <<
G4endl;
230 while( (i<9) && (ran>=kkb[i]) ) {
231 ++i;
232 loop++;
233 if (loop > 1000) {
235 break;
236 }
237 }
238
239 if( i == 9 )return true;
240
241
242
243
244 switch( ipakkb1[i] )
245 {
246 case 10:
247 vec[i3]->SetDefinition( aKaonPlus );
248 vec[i3]->SetMayBeKilled(false);
249 break;
250 case 11:
251 vec[i3]->SetDefinition( aKaonZS );
252 vec[i3]->SetMayBeKilled(false);
253 break;
254 case 12:
255 vec[i3]->SetDefinition( aKaonZL );
256 vec[i3]->SetMayBeKilled(false);
257 break;
258 }
259
260 if( vecLen == 1 )
261 {
263 switch( ipakkb2[i] )
264 {
265 case 11:
268 break;
269 case 12:
272 break;
273 case 13:
276 break;
277 }
280
281 }
282 else
283 {
284 switch( ipakkb2[i] )
285 {
286 case 11:
287 vec[i4]->SetDefinition( aKaonZS );
288 vec[i4]->SetMayBeKilled(false);
289 break;
290 case 12:
291 vec[i4]->SetDefinition( aKaonZL );
292 vec[i4]->SetMayBeKilled(false);
293 break;
294 case 13:
295 vec[i4]->SetDefinition( aKaonMinus );
296 vec[i4]->SetMayBeKilled(false);
297 break;
298 }
299 }
300
301 } else if( ran < avy ) {
302 if( availableEnergy < 1.6 )return true;
303
304 const G4double ky[] = { 0.200, 0.300, 0.400, 0.550, 0.625, 0.700,
305 0.800, 0.850, 0.900, 0.950, 0.975, 1.000 };
306 const G4int ipaky1[] = { 18, 18, 18, 20, 20, 20, 21, 21, 21, 22, 22, 22 };
307 const G4int ipaky2[] = { 10, 11, 12, 10, 11, 12, 10, 11, 12, 10, 11, 12 };
308 const G4int ipakyb1[] = { 19, 19, 19, 23, 23, 23, 24, 24, 24, 25, 25, 25 };
309 const G4int ipakyb2[] = { 13, 12, 11, 13, 12, 11, 13, 12, 11, 13, 12, 11 };
311 i = 0;
312
313 loop = 0;
315 edb <<
" While count exceeded " <<
G4endl;
316 while( (i<12) && (ran>ky[i]) ) {
317 ++i;
318 loop++;
319 if (loop > 1000) {
321 break;
322 }
323 }
324
325 if( i == 12 )return true;
327
328
329
330
331
332 switch( ipaky1[i] )
333 {
334 case 18:
336 break;
337 case 20:
339 break;
340 case 21:
342 break;
343 case 22:
345 break;
346 }
347 targetHasChanged = true;
348 switch( ipaky2[i] )
349 {
350 case 10:
351 vec[i3]->SetDefinition( aKaonPlus );
352 vec[i3]->SetMayBeKilled(false);
353 break;
354 case 11:
355 vec[i3]->SetDefinition( aKaonZS );
356 vec[i3]->SetMayBeKilled(false);
357 break;
358 case 12:
359 vec[i3]->SetDefinition( aKaonZL );
360 vec[i3]->SetMayBeKilled(false);
361 break;
362 }
363
364 } else {
365
366
370 (currentMass > sigmaMinusMass) ) {
371 switch( ipakyb1[i] )
372 {
373 case 19:
375 break;
376 case 23:
378 break;
379 case 24:
381 break;
382 case 25:
384 break;
385 }
386 incidentHasChanged = true;
387 switch( ipakyb2[i] )
388 {
389 case 11:
390 vec[i3]->SetDefinition( aKaonZS );
391 vec[i3]->SetMayBeKilled(false);
392 break;
393 case 12:
394 vec[i3]->SetDefinition( aKaonZL );
395 vec[i3]->SetMayBeKilled(false);
396 break;
397 case 13:
398 vec[i3]->SetDefinition( aKaonMinus );
399 vec[i3]->SetMayBeKilled(false);
400 break;
401 }
402
403 } else {
404 switch( ipaky1[i] )
405 {
406 case 18:
408 break;
409 case 20:
411 break;
412 case 21:
414 break;
415 case 22:
417 break;
418 }
419 incidentHasChanged = true;
420 switch( ipaky2[i] )
421 {
422 case 10:
423 vec[i3]->SetDefinition( aKaonPlus );
424 vec[i3]->SetMayBeKilled(false);
425 break;
426 case 11:
427 vec[i3]->SetDefinition( aKaonZS );
428 vec[i3]->SetMayBeKilled(false);
429 break;
430 case 12:
431 vec[i3]->SetDefinition( aKaonZL );
432 vec[i3]->SetMayBeKilled(false);
433 break;
434 }
435 }
436 }
437 }
438 else return true;
439
440
441
442
443
444
445
446
447
448
449 currentMass = currentParticle.
GetMass()/GeV;
450 targetMass = targetParticle.
GetMass()/GeV;
451
452 G4double energyCheck = centerofmassEnergy-(currentMass+targetMass);
453 for( i=0; i<vecLen; ++i )
454 {
455 energyCheck -= vec[i]->GetMass()/GeV;
456 if( energyCheck < 0.0 )
457 {
458 vecLen = std::max( 0, --i );
460 for(j=i; j<vecLen; j++) delete vec[j];
461 break;
462 }
463 }
464
465 return true;
466}
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::ostringstream G4ExceptionDescription
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4double G4Log(G4double x)
static G4AntiLambda * AntiLambda()
static G4AntiNeutron * AntiNeutron()
static G4AntiProton * AntiProton()
static G4AntiSigmaMinus * AntiSigmaMinus()
static G4AntiSigmaPlus * AntiSigmaPlus()
static G4AntiSigmaZero * AntiSigmaZero()
G4ParticleDefinition * GetDefinition() const
void SetElement(G4int anIndex, Type *anElement)
static G4KaonMinus * KaonMinus()
static G4KaonPlus * KaonPlus()
static G4KaonZeroLong * KaonZeroLong()
static G4KaonZeroShort * KaonZeroShort()
static G4Lambda * Lambda()
static G4Neutron * Neutron()
G4double GetPDGMass() const
static G4Proton * Proton()
const G4ParticleDefinition * GetDefinition() const
G4double GetTotalEnergy() const
void SetMayBeKilled(const G4bool f)
void SetSide(const G4int sid)
void SetDefinitionAndUpdateE(const G4ParticleDefinition *aParticleDefinition)
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
static G4SigmaMinus * SigmaMinus()
static G4SigmaPlus * SigmaPlus()
static G4SigmaZero * SigmaZero()