226{
227
228
229
230
231
232
235
236
237
238
239
250
251
252
253
260 {
261 G4cout <<
"########################################"
262 <<"########################################"
266 G4cout <<
"Initial projectile A=" <<AP
267 <<", Z=" <<ZP
268 <<", radius = " <<rP/fermi <<" fm"
270 G4cout <<
"Initial target A=" <<AT
271 <<", Z=" <<ZT
272 <<", radius = " <<rT/fermi <<" fm"
274 G4cout <<
"Projectile momentum and Energy/nuc = " <<pP <<
" ," <<E <<
G4endl;
275 }
276
277
278
279
280
281 G4double rm = ZP * ZT * elm_coupling / (E * AP);
284
285
286
287
288
289
294
295
296
297
298
299 G4bool skipInteraction =
false;
300 const G4int maxNumberOfLoops = 1000;
301 G4int loopCounter = -1;
302 while (Dabr == 0 && ++loopCounter < maxNumberOfLoops)
303 {
304
305
306
307
308
309
312
313
314
315
316
317
318
319 if (rm >= fradius * rPT) {
320 skipInteraction = true;
321 }
322
323
324
325
326
328 r = 1.1 * rPT;
329 while (r > rPT && ++evtcnt < 1000)
330 {
332 r = (rm + std::sqrt(rm*rm + 4.0*bsq)) / 2.0;
333 }
334
335
336
337
338 if (evtcnt >= 1000) {
339 skipInteraction = true;
340 }
341
342 rsq = r * r;
343
344
345
346
347 if (rT > rP)
348 {
349 G4double x = (rPsq + rsq - rTsq) / 2.0 / r;
350 if (x > 0.0) CT = 2.0 * std::sqrt(rTsq - x*x);
351 else CT = 2.0 * std::sqrt(rTsq - rsq);
352 }
353 else
354 {
355 G4double x = (rTsq + rsq - rPsq) / 2.0 / r;
356 if (x > 0.0) CT = 2.0 * std::sqrt(rTsq - x*x);
357 else CT = 2.0 * rT;
358 }
359
360
361
362
363
364
365
366
367 delete theAbrasionGeometry;
369 F = theAbrasionGeometry->
F();
373 for (
G4int i = 0; i<10; ++i)
374 {
376 if (n > 0)
377 {
378 if (n>AP) Dabr = (
G4int) AP;
379 else Dabr = (
G4int) n;
380 break;
381 }
382 }
383 }
384
385 if ( loopCounter >= maxNumberOfLoops || skipInteraction ) {
386
391 G4cout <<
"Particle energy too low to overcome repulsion." <<
G4endl;
392 G4cout <<
"Event rejected and original track maintained" <<
G4endl;
393 G4cout <<
"########################################"
394 <<"########################################"
396 }
397 delete theAbrasionGeometry;
399 }
400
402 {
404 G4cout <<
"Impact parameter = " <<r/fermi <<
" fm" <<
G4endl;
406 }
407
408
409
410
411
412
413
414
417
418
419
420
421
422
423
425 G4Fragment *fragmentP = GetAbradedNucleons (Dabr, AP, ZP, rP);
428 for (i=0; i<nSecP; ++i)
429 {
431 GetParticle()->GetTotalEnergy();
432 }
433
434
435
436
437
439 if (DspcP <= 0) DspcP = 0;
440 else if (DspcP > AP-Dabr) DspcP = ((
G4int) AP) - Dabr;
441
442
443
444
445
446
447
448 G4bool excitationAbsorbedByProjectile =
false;
449 if (fragmentP != nullptr)
450 {
453 if (Dabr < AT)
455 if (excitationAbsorbedByProjectile)
456 ExP = GetNucleonInducedExcitation(rP, rT, r);
458 if (xP > B*(AP-Dabr)) xP = B*(AP-Dabr);
460 lorentzVector.
setE(lorentzVector.e()+xP);
462 TotalEPost += lorentzVector.e();
463 }
465
466
467
468
469
470
471
472
473
474 G4Fragment *fragmentT = GetAbradedNucleons (Dabr, AT, ZT, rT);
476 for (i=nSecP; i<nSec; ++i)
477 {
479 GetParticle()->GetTotalEnergy();
480 }
481
482
483
484
485
487 if (DspcT <= 0) DspcT = 0;
488 else if (DspcT > AP-Dabr) DspcT = ((
G4int) AT) - Dabr;
489
490
491
492
493
494
495
496 if (fragmentT != nullptr)
497 {
500 if (!excitationAbsorbedByProjectile)
501 ExT = GetNucleonInducedExcitation(rT, rP, r);
503 if (xT > B*(AT-Dabr)) xT = B*(AT-Dabr);
505 lorentzVector.
setE(lorentzVector.e()+xT);
507 TotalEPost += lorentzVector.e();
508 }
509
510
511
512
513
514
515 G4double deltaE = TotalEPre - TotalEPost;
516 if (deltaE > 0.0 && conserveEnergy)
517 {
519 boost = boost / boost.
mag() * beta;
520 }
521
522
523
524
526 for (i=0; i<nSecP; ++i)
527 {
529 GetParticle();
531 lorentzVector.
boost(-boost);
533 pBalance -= lorentzVector.
vect();
534 }
535
536
537
538
539
540
541
542
543
544
545 if (fragmentP != nullptr)
546 {
549 if (conserveMomentum)
552 else
553 {
555 fragmentP->
SetMomentum(lorentzVector.
boost(-boost * fragmentGroundStateM/fragmentM));
556 }
557 }
558
559
560
561
563 {
565 G4cout <<
"-----------------------------------" <<
G4endl;
566 G4cout <<
"Secondary nucleons from projectile:" <<
G4endl;
567 G4cout <<
"-----------------------------------" <<
G4endl;
569 for (i=0; i<nSecP; ++i)
570 {
577 }
581 if (fragmentP != nullptr)
583 else
590 for (i=nSecP; i<nSec; ++i)
591 {
598 }
602 if (fragmentT != nullptr)
604 else
606 }
607
608
609
610
611
612 if (fragmentP !=nullptr)
613 {
615
616 products = theExcitationHandler->
BreakItUp(*fragmentP);
617
618
619 delete fragmentP;
620 fragmentP = nullptr;
621
622 G4ReactionProductVector::iterator iter;
623 for (iter = products->begin(); iter != products->end(); ++iter)
624 {
627 (*iter)->GetTotalEnergy(), (*iter)->GetMomentum());
629 G4String particleName = (*iter)->GetDefinition()->GetParticleName();
630 delete (*iter);
631 if (
verboseLevel >= 2 && particleName.find(
"[",0) < particleName.size())
632 {
636 G4cout <<
" fragmentP = " <<particleName
639 }
640 }
641 delete products;
642 }
643
644
645
646
647
648
649 if (fragmentT != nullptr)
650 {
652
653 products = theExcitationHandler->
BreakItUp(*fragmentT);
654
655
656
657 fragmentT = nullptr;
658
659 G4ReactionProductVector::iterator iter;
660 for (iter = products->begin(); iter != products->end(); ++iter)
661 {
664 (*iter)->GetTotalEnergy(), (*iter)->GetMomentum());
666 G4String particleName = (*iter)->GetDefinition()->GetParticleName();
667 delete (*iter);
668 if (
verboseLevel >= 2 && particleName.find(
"[",0) < particleName.size())
669 {
673 G4cout <<
" fragmentT = " <<particleName
676 }
677 }
678 delete products;
679 }
680
682 G4cout <<
"########################################"
683 <<"########################################"
685
686 delete theAbrasionGeometry;
688}
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 AtomicMass(const G4double A, const G4double Z, const G4int numberOfLambdas=0) const
G4double GetEnergyDeposit()
G4double GetPDGCharge() const
G4int GetBaryonNumber() const
const G4String & GetParticleName() const
static G4Pow * GetInstance()
G4double powA(G4double A, G4double y) const
G4double GetWilsonRadius(G4double A)