217{
218
219
220
221
222
223
226
227
228
229
230
241
242
243
244
251 {
252 G4cout <<
"########################################"
253 <<"########################################"
257 G4cout <<
"Initial projectile A=" <<AP
258 <<", Z=" <<ZP
259 <<", radius = " <<rP/fermi <<" fm"
261 G4cout <<
"Initial target A=" <<AT
262 <<", Z=" <<ZT
263 <<", radius = " <<rT/fermi <<" fm"
265 G4cout <<
"Projectile momentum and Energy/nuc = " <<pP <<
" ," <<E <<
G4endl;
266 }
267
268
269
270
271
272 G4double rm = ZP * ZT * elm_coupling / (E * AP);
275
276
277
278
279
280
285
286
287
288
289
290 G4bool skipInteraction =
false;
291 const G4int maxNumberOfLoops = 1000;
292 G4int loopCounter = -1;
293 while (Dabr == 0 && ++loopCounter < maxNumberOfLoops)
294 {
295
296
297
298
299
300
303
304
305
306
307
308
309
310 if (rm >= fradius * rPT) {
311 skipInteraction = true;
312 }
313
314
315
316
317
319 r = 1.1 * rPT;
320 while (r > rPT && ++evtcnt < 1000)
321 {
323 r = (rm + std::sqrt(rm*rm + 4.0*bsq)) / 2.0;
324 }
325
326
327
328
329 if (evtcnt >= 1000) {
330 skipInteraction = true;
331 }
332
333 rsq = r * r;
334
335
336
337
338 if (rT > rP)
339 {
340 G4double x = (rPsq + rsq - rTsq) / 2.0 / r;
341 if (x > 0.0) CT = 2.0 * std::sqrt(rTsq - x*x);
342 else CT = 2.0 * std::sqrt(rTsq - rsq);
343 }
344 else
345 {
346 G4double x = (rTsq + rsq - rPsq) / 2.0 / r;
347 if (x > 0.0) CT = 2.0 * std::sqrt(rTsq - x*x);
348 else CT = 2.0 * rT;
349 }
350
351
352
353
354
355
356
357
358 delete theAbrasionGeometry;
360 F = theAbrasionGeometry->
F();
364 for (
G4int i = 0; i<10; ++i)
365 {
367 if (n > 0)
368 {
369 if (n>AP) Dabr = (
G4int) AP;
370 else Dabr = (
G4int) n;
371 break;
372 }
373 }
374 }
375
376 if ( loopCounter >= maxNumberOfLoops || skipInteraction ) {
377
382 G4cout <<
"Particle energy too low to overcome repulsion." <<
G4endl;
383 G4cout <<
"Event rejected and original track maintained" <<
G4endl;
384 G4cout <<
"########################################"
385 <<"########################################"
387 }
388 delete theAbrasionGeometry;
390 }
391
393 {
395 G4cout <<
"Impact parameter = " <<r/fermi <<
" fm" <<
G4endl;
397 }
398
399
400
401
402
403
404
405
408
409
410
411
412
413
414
416 G4Fragment *fragmentP = GetAbradedNucleons (Dabr, AP, ZP, rP);
419 for (i=0; i<nSecP; ++i)
420 {
422 GetParticle()->GetTotalEnergy();
423 }
424
425
426
427
428
430 if (DspcP <= 0) DspcP = 0;
431 else if (DspcP > AP-Dabr) DspcP = ((
G4int) AP) - Dabr;
432
433
434
435
436
437
438
439 G4bool excitationAbsorbedByProjectile =
false;
440 if (fragmentP != nullptr)
441 {
444 if (Dabr < AT)
446 if (excitationAbsorbedByProjectile)
447 ExP = GetNucleonInducedExcitation(rP, rT, r);
449 if (xP > B*(AP-Dabr)) xP = B*(AP-Dabr);
451 lorentzVector.
setE(lorentzVector.e()+xP);
453 TotalEPost += lorentzVector.e();
454 }
456
457
458
459
460
461
462
463
464
465 G4Fragment *fragmentT = GetAbradedNucleons (Dabr, AT, ZT, rT);
467 for (i=nSecP; i<nSec; ++i)
468 {
470 GetParticle()->GetTotalEnergy();
471 }
472
473
474
475
476
478 if (DspcT <= 0) DspcT = 0;
479 else if (DspcT > AP-Dabr) DspcT = ((
G4int) AT) - Dabr;
480
481
482
483
484
485
486
487 if (fragmentT != nullptr)
488 {
491 if (!excitationAbsorbedByProjectile)
492 ExT = GetNucleonInducedExcitation(rT, rP, r);
494 if (xT > B*(AT-Dabr)) xT = B*(AT-Dabr);
496 lorentzVector.
setE(lorentzVector.e()+xT);
498 TotalEPost += lorentzVector.e();
499 }
500
501
502
503
504
505
506 G4double deltaE = TotalEPre - TotalEPost;
507 if (deltaE > 0.0 && conserveEnergy)
508 {
510 boost = boost / boost.
mag() * beta;
511 }
512
513
514
515
517 for (i=0; i<nSecP; ++i)
518 {
520 GetParticle();
522 lorentzVector.
boost(-boost);
524 pBalance -= lorentzVector.
vect();
525 }
526
527
528
529
530
531
532
533
534
535
536 if (fragmentP != nullptr)
537 {
540 if (conserveMomentum)
543 else
544 {
546 fragmentP->
SetMomentum(lorentzVector.
boost(-boost * fragmentGroundStateM/fragmentM));
547 }
548 }
549
550
551
552
554 {
556 G4cout <<
"-----------------------------------" <<
G4endl;
557 G4cout <<
"Secondary nucleons from projectile:" <<
G4endl;
558 G4cout <<
"-----------------------------------" <<
G4endl;
560 for (i=0; i<nSecP; ++i)
561 {
568 }
572 if (fragmentP != nullptr)
574 else
581 for (i=nSecP; i<nSec; ++i)
582 {
589 }
593 if (fragmentT != nullptr)
595 else
597 }
598
599
600
601
602
603 if (fragmentP !=nullptr)
604 {
606
607 products = theExcitationHandler->
BreakItUp(*fragmentP);
608
609
610 delete fragmentP;
611 fragmentP = nullptr;
612
613 G4ReactionProductVector::iterator iter;
614 for (iter = products->begin(); iter != products->end(); ++iter)
615 {
618 (*iter)->GetTotalEnergy(), (*iter)->GetMomentum());
620 G4String particleName = (*iter)->GetDefinition()->GetParticleName();
621 delete (*iter);
622 if (
verboseLevel >= 2 && particleName.find(
"[",0) < particleName.size())
623 {
627 G4cout <<
" fragmentP = " <<particleName
630 }
631 }
632 delete products;
633 }
634
635
636
637
638
639
640 if (fragmentT != nullptr)
641 {
643
644 products = theExcitationHandler->
BreakItUp(*fragmentT);
645
646
647
648 fragmentT = nullptr;
649
650 G4ReactionProductVector::iterator iter;
651 for (iter = products->begin(); iter != products->end(); ++iter)
652 {
655 (*iter)->GetTotalEnergy(), (*iter)->GetMomentum());
657 G4String particleName = (*iter)->GetDefinition()->GetParticleName();
658 delete (*iter);
659 if (
verboseLevel >= 2 && particleName.find(
"[",0) < particleName.size())
660 {
664 G4cout <<
" fragmentT = " <<particleName
667 }
668 }
669 delete products;
670 }
671
673 G4cout <<
"########################################"
674 <<"########################################"
676
677 delete theAbrasionGeometry;
679}
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
CLHEP::HepLorentzVector G4LorentzVector
G4long G4Poisson(G4double mean)
std::vector< G4ReactionProduct * > G4ReactionProductVector
G4GLOB_DLL std::ostream G4cout
HepLorentzVector & boost(double, double, double)
Hep3Vector findBoostToCM() const
void DumpInfo(G4int mode=0) const
G4ParticleDefinition * GetDefinition() const
G4LorentzVector Get4Momentum() const
G4double GetKineticEnergy() const
void Set4Momentum(const G4LorentzVector &momentum)
G4ReactionProductVector * BreakItUp(const G4Fragment &theInitialState)
G4double GetGroundStateMass() const
const G4LorentzVector & GetMomentum() const
void SetMomentum(const G4LorentzVector &value)
void SetStatusChange(G4HadFinalStateStatus aS)
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
std::size_t GetNumberOfSecondaries() const
void SetEnergyChange(G4double anEnergy)
G4HadSecondary * GetSecondary(size_t i)
void SetMomentumChange(const G4ThreeVector &aV)
const G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
const G4LorentzVector & Get4Momentum() const
G4double GetTotalEnergy() const
G4DynamicParticle * GetParticle()
G4HadFinalState theParticleChange
G4double GetExcitationEnergyOfTarget()
G4double GetExcitationEnergyOfProjectile()
G4double GetEnergyDeposit()
G4double AtomicMass(const G4double A, const G4double Z) const
G4double GetPDGCharge() const
G4int GetBaryonNumber() const
const G4String & GetParticleName() const
static G4Pow * GetInstance()
G4double powA(G4double A, G4double y) const
G4double GetWilsonRadius(G4double A)