243{
244
245
246
247
248
249
252
253
254
255
256
267
268
269
270
277 {
278 G4cout <<
"########################################"
279 <<"########################################"
283 G4cout <<
"Initial projectile A=" <<AP
284 <<", Z=" <<ZP
285 <<", radius = " <<rP/fermi <<" fm"
287 G4cout <<
"Initial target A=" <<AT
288 <<", Z=" <<ZT
289 <<", radius = " <<rT/fermi <<" fm"
291 G4cout <<
"Projectile momentum and Energy/nuc = " <<pP <<
" ," <<E <<
G4endl;
292 }
293
294
295
296
297
298 G4double rm = ZP * ZT * elm_coupling / (E * AP);
301
302
303
304
305
306
311
312
313
314
315
316 while (Dabr == 0)
317 {
318
319 if (theAbrasionGeometry)
320 {
321 delete theAbrasionGeometry;
322 theAbrasionGeometry = NULL;
323 }
324
325
326
327
328
329
332
333
334
335
336
337
338
339
340 if (rm >= fradius * rPT) {
345 G4cout <<
"Particle energy too low to overcome repulsion." <<
G4endl;
346 G4cout <<
"Event rejected and original track maintained" <<
G4endl;
347 G4cout <<
"########################################"
348 <<"########################################"
350 }
352 }
353
354
355
356
357
359 r = 1.1 * rPT;
360 while (r > rPT && ++evtcnt < 1000)
361 {
363 r = (rm + std::sqrt(rm*rm + 4.0*bsq)) / 2.0;
364 }
365
366
367
368
369
370 if (evtcnt >= 1000) {
375 G4cout <<
"Particle energy too low to overcome repulsion." <<
G4endl;
376 G4cout <<
"Event rejected and original track maintained" <<
G4endl;
377 G4cout <<
"########################################"
378 <<"########################################"
380 }
382 }
383
384
385 rsq = r * r;
386
387
388
389
390 if (rT > rP)
391 {
392 G4double x = (rPsq + rsq - rTsq) / 2.0 / r;
393 if (x > 0.0) CT = 2.0 * std::sqrt(rTsq - x*x);
394 else CT = 2.0 * std::sqrt(rTsq - rsq);
395 }
396 else
397 {
398 G4double x = (rTsq + rsq - rPsq) / 2.0 / r;
399 if (x > 0.0) CT = 2.0 * std::sqrt(rTsq - x*x);
400 else CT = 2.0 * rT;
401 }
402
403
404
405
406
407
408
409
411 F = theAbrasionGeometry->
F();
413 G4double Mabr = F * AP * (1.0 - std::exp(-CT/lambda));
415 for (
G4int i = 0; i<10; i++)
416 {
418 if (n > 0)
419 {
420 if (n>AP) Dabr = (
G4int) AP;
421 else Dabr = (
G4int) n;
422 break;
423 }
424 }
425 }
427 {
429 G4cout <<
"Impact parameter = " <<r/fermi <<
" fm" <<
G4endl;
431 }
432
433
434
435
436
437
438
439
442
443
444
445
446
447
448
450 G4Fragment *fragmentP = GetAbradedNucleons (Dabr, AP, ZP, rP);
453 for (i=0; i<nSecP; i++)
454 {
456 GetParticle()->GetTotalEnergy();
457 }
458
459
460
461
462
464 if (DspcP <= 0) DspcP = 0;
465 else if (DspcP > AP-Dabr) DspcP = ((
G4int) AP) - Dabr;
466
467
468
469
470
471
472
473 G4bool excitationAbsorbedByProjectile =
false;
474 if (fragmentP != NULL)
475 {
478 if (Dabr < AT)
480 if (excitationAbsorbedByProjectile)
481 ExP = GetNucleonInducedExcitation(rP, rT, r);
483 if (xP > B*(AP-Dabr)) xP = B*(AP-Dabr);
485 lorentzVector.
setE(lorentzVector.e()+xP);
487 TotalEPost += lorentzVector.e();
488 }
490
491
492
493
494
495
496
497
498
499 G4Fragment *fragmentT = GetAbradedNucleons (Dabr, AT, ZT, rT);
501 for (i=nSecP; i<nSec; i++)
502 {
504 GetParticle()->GetTotalEnergy();
505 }
506
507
508
509
510
512 if (DspcT <= 0) DspcT = 0;
513 else if (DspcT > AP-Dabr) DspcT = ((
G4int) AT) - Dabr;
514
515
516
517
518
519
520
521 if (fragmentT != NULL)
522 {
525 if (!excitationAbsorbedByProjectile)
526 ExT = GetNucleonInducedExcitation(rT, rP, r);
528 if (xT > B*(AT-Dabr)) xT = B*(AT-Dabr);
530 lorentzVector.
setE(lorentzVector.e()+xT);
532 TotalEPost += lorentzVector.e();
533 }
534
535
536
537
538
539
540 G4double deltaE = TotalEPre - TotalEPost;
541 if (deltaE > 0.0 && conserveEnergy)
542 {
543 G4double beta = std::sqrt(1.0 - EMassP*EMassP/std::pow(deltaE+EMassP,2.0));
544 boost = boost / boost.
mag() * beta;
545 }
546
547
548
549
551 for (i=0; i<nSecP; i++)
552 {
554 GetParticle();
556 lorentzVector.
boost(-boost);
558 pBalance -= lorentzVector.
vect();
559 }
560
561
562
563
564
565
566
567
568
569
570 if (fragmentP != NULL)
571 {
574 if (conserveMomentum)
577 else
578 {
580 fragmentP->
SetMomentum(lorentzVector.
boost(-boost * fragmentGroundStateM/fragmentM));
581 }
582 }
583
584
585
586
588 {
590 G4cout <<
"-----------------------------------" <<
G4endl;
591 G4cout <<
"Secondary nucleons from projectile:" <<
G4endl;
592 G4cout <<
"-----------------------------------" <<
G4endl;
594 for (i=0; i<nSecP; i++)
595 {
602 }
606 if (fragmentP != NULL)
608 else
615 for (i=nSecP; i<nSec; i++)
616 {
623 }
627 if (fragmentT != NULL)
629 else
631 }
632
633
634
635
636
637 if (fragmentP !=NULL)
638 {
640 if (fragmentP->
GetZ() != fragmentP->
GetA())
641 products = theExcitationHandler->
BreakItUp(*fragmentP);
642 else
643 products = theExcitationHandlerx->
BreakItUp(*fragmentP);
644 delete fragmentP;
645 fragmentP = NULL;
646
647 G4ReactionProductVector::iterator iter;
648 for (iter = products->begin(); iter != products->end(); ++iter)
649 {
652 (*iter)->GetTotalEnergy(), (*iter)->GetMomentum());
654 G4String particleName = (*iter)->GetDefinition()->GetParticleName();
655 delete (*iter);
656 if (
verboseLevel >= 2 && particleName.find(
"[",0) < particleName.size())
657 {
661 G4cout <<
" fragmentP = " <<particleName
664 }
665 }
666 delete products;
667 }
668
669
670
671
672
673
674 if (fragmentT != NULL)
675 {
677 if (fragmentT->
GetZ() != fragmentT->
GetA())
678 products = theExcitationHandler->
BreakItUp(*fragmentT);
679 else
680 products = theExcitationHandlerx->
BreakItUp(*fragmentT);
681 delete fragmentT;
682 fragmentT = NULL;
683
684 G4ReactionProductVector::iterator iter;
685 for (iter = products->begin(); iter != products->end(); ++iter)
686 {
689 (*iter)->GetTotalEnergy(), (*iter)->GetMomentum());
691 G4String particleName = (*iter)->GetDefinition()->GetParticleName();
692 delete (*iter);
693 if (
verboseLevel >= 2 && particleName.find(
"[",0) < particleName.size())
694 {
698 G4cout <<
" fragmentT = " <<particleName
701 }
702 }
703 delete products;
704 }
705
707 G4cout <<
"########################################"
708 <<"########################################"
710
711 delete theAbrasionGeometry;
712
714}
CLHEP::HepLorentzVector G4LorentzVector
G4long G4Poisson(G4double mean)
std::vector< G4ReactionProduct * > G4ReactionProductVector
G4DLLIMPORT 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) const
G4double GetGroundStateMass() const
const G4LorentzVector & GetMomentum() const
void SetMomentum(const G4LorentzVector &value)
void SetStatusChange(G4HadFinalStateStatus aS)
G4int GetNumberOfSecondaries() const
void AddSecondary(G4DynamicParticle *aP)
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
G4double GetWilsonRadius(G4double A)