531{
532 G4double eMass = CLHEP::electron_mass_c2;
534
535
541
545
546
548
549 G4double mingEnergy = std::sqrt(Q2);
551 G4double intervalgEnergy = maxgEnergy - mingEnergy;
552
553
555
556
557 G4double gMomentum = std::sqrt(gEnergy*gEnergy - Q2);
558
559
560 G4double particleFinalEnergy = kinEnergy - gEnergy;
561 G4double particleFinalMomentum = std::sqrt(particleFinalEnergy*particleFinalEnergy -
563
568
569
570
571
572
573
575
576
579 G4double maxEnergy = std::min(tmax, maxPairEnergy);
581
582 if (minEnergy >= maxEnergy) { return; }
583
584
585
586
587
588 CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
589
591
592
595
596
597
598
600
601
602
603
604 G4int iz1(0), iz2(0);
607 iz1 = iz2 = iz;
608 break;
610 iz2 = iz;
611 if(iz > 0) { iz1 = iz-1; }
612 else { iz1 = iz2; }
613 break;
614 }
615 }
616 if (0 == iz1) { iz1 = iz2 =
NZDATPAIR-1; }
617
620
621 do {
622 ++count;
623
625
627 if(iz1 != iz2) {
631
632
633 x += (x2 - x)*(
lnZ - lz1)/(lz2 - lz1);
634 }
635
636 pairEnergy = kinEnergy*
G4Exp(x*coeff);
637
638
639 } while((pairEnergy < minEnergy || pairEnergy > maxEnergy) && 50 > count);
640
641
642
645 fAngularGenerator->PhiRotation(partDirection, phi3);
646
650
652 G4double A1 = -(Q2 - 2.*kinEnergy*gEnergy);
653 G4double B1 = -(2.*gMomentum*particleMomentum);
655
657 G4double sintg = std::sqrt((1.0 - costg)*(1.0 + costg));
661
663 dirGamma.
set(sintg*cospg, sintg*sinpg, costg);
665
666 G4double Ap = particleMomentum*particleMomentum +
667 particleFinalMomentum*particleFinalMomentum + gMomentum*gMomentum;
668 G4double A = Ap - 2.*particleMomentum*gMomentum*costg;
669 G4double B = 2.*particleMomentum*gMomentum*sintg*cospg;
670 G4double C = 2.*particleFinalMomentum*gMomentum*costg -
671 2.*particleMomentum*particleFinalMomentum;
673 G4double t1interval = (1./(
A +
C + absB*mint3) - 1./(
A +
C + absB*maxt3))/absB;
677
678
680 partDirection.
set(sint1, 0., cost1);
681 fAngularGenerator->PhiRotation(partDirection, Phi);
682 kinEnergy = particleFinalEnergy;
683
686
687 G4LorentzVector eFourMomentumMQ = fAngularGenerator->eDP2(Q2, eMass2, eMass2, cost5, phi5);
688 G4LorentzVector pFourMomentumMQ = fAngularGenerator->pDP2(eMass2, eFourMomentumMQ);
689
692
693 eEnergy = eFourMomentum.
t();
694 pEnergy = pFourMomentum.
t();
695
697 G4double A3 = Q2 + 2.*gEnergy*particleFinalEnergy;
698 G4double B3 = -2.*gMomentum*particleFinalMomentum;
699
703
704 G4double sintQ3 = std::sqrt(1. - tQMG*tQMG);
707
708 G4double Ap = particleMomentum*particleMomentum +
709 particleFinalMomentum*particleFinalMomentum + gMomentum*gMomentum;
710 G4double A = Ap + 2.*particleFinalMomentum*gMomentum*tQMG;
711 G4double B = -2.*particleMomentum*gMomentum*sintQ3*cospQP;
712 G4double C = -2.*particleMomentum*gMomentum*tQMG - 2.*particleMomentum*particleFinalMomentum;
713
715 G4double t3interval = (1./(
A +
C + absB*mint3) - 1./(
A +
C + absB*maxt3))/absB;
719
720 G4double cost = -sint3*sintQ3*cospQP + cost3*tQMG;
721 G4double sint = std::sqrt((1. + cost)*(1. - cost));
722 G4double cosp = (sintQ3*cospQP*cost3 + sint3*tQMG)/sint;
724
726 dirGamma.
set(sint*cosp, sint*sinp, cost);
728
729
731 partDirection.
set(sint3, 0., cost3);
732 fAngularGenerator->PhiRotation(partDirection, Phi);
733 kinEnergy = particleFinalEnergy;
734
737
738 G4LorentzVector eFourMomentumMQ = fAngularGenerator->eDP2(Q2, eMass2, eMass2, cost5, phi5);
739 G4LorentzVector pFourMomentumMQ = fAngularGenerator->pDP2(eMass2, eFourMomentumMQ);
740
743
744 eEnergy = eFourMomentum.
t();
745 pEnergy = pFourMomentum.
t();
746
752 particleFinalMomentum = std::sqrt(particleFinalEnergy*particleFinalEnergy -
754
756 G4double maxeEnergy = kinEnergy - particleFinalEnergy - eMass;
757 G4double eEnergyinterval = maxeEnergy - mineEnergy;
758 eEnergy = mineEnergy + eEnergyinterval*
randNumbs[4];
759
766
767 G4double eMomentum = std::sqrt(eEnergy*eEnergy - eMass*eMass);
768 pEnergy = kinEnergy - particleFinalEnergy - eEnergy;
769 G4double pMomentum = std::sqrt(pEnergy*pEnergy - eMass*eMass);
770
772 G4double B3 = -2.*particleMomentum*particleFinalMomentum;
773 G4double cost3interval =
G4Log((A3 + B3*Cmax)/(A3 + B3*Cmin))/B3;
775 G4double cost3 = A3*(expanCost3r6 - 1.)/B3 + Cmin*expanCost3r6;
776 G4double sint3 = std::sqrt((1. - cost3)*(1. + cost3));
777
778 partDirection.
set(sint3, 0., cost3);
779
781 muFinalMomentumVector.
set(particleFinalMomentum*sint3, 0., particleFinalMomentum*cost3);
782
784 G4LorentzVector muFinalFourMomentum(muFinalMomentumVector, particleFinalEnergy);
786 G4double A5 = auxVec1.
mag2() - 2.*eEnergy*(kinEnergy - particleFinalEnergy) +
787 2.*particleMomentumVector[2]*eMomentum - 2.*particleFinalMomentum*eMomentum*cost3;
788 G4double B5 = -2.*particleFinalMomentum*eMomentum*(sint3*cosp3*cosp5 + sint3*sinp3*sinp5);
793 G4double t5interval =
G4Log((absA5 + absB5*maxt5)/(absA5 + absB5*mint5))/absB5;
798
799 eDirection.
set(sint5*cosp5, sint5*sinp5, cost5);
801
802 G4ThreeVector auxVec2 = particleMomentumVector - muFinalMomentumVector - eMomentumVector;
804 G4double Bp = particleFinalMomentum*(sint3*cosp3*cosp6 + sint3*sinp3*sinp6) +
805 eMomentum*(sint5*cosp5*cosp6 + sint5*sinp5*sinp6);
806 G4double Cp = -particleMomentum + particleFinalMomentum*cost3 + eMomentum*cost5;
807 G4double A6 = p1mp3mp52 + pMomentum*pMomentum;
812 G4double absA6C6 = std::abs(A6 + C6);
814 G4double t6interval = (1./(absA6C6 + absB6*mint6) - 1./(absA6C6 + absB6*maxt6))/absB6;
815 G4double t6 = (-absA6C6 + 1./(1./(absA6C6 + absB6*mint6) - absB6*t6interval*
randNumbs[8]))/absB6;
818
819 pDirection.
set(sint6*cosp6, sint6*sinp6, cost6);
820
821 } else {
826 particleFinalMomentum = std::sqrt(particleFinalEnergy*particleFinalEnergy -
828
829 G4double maxpEnergy = kinEnergy - particleFinalEnergy - eMass;
830 G4double pEnergyinterval = maxpEnergy - eMass;
831 pEnergy = eMass + pEnergyinterval*
randNumbs[4];
832
839
840 G4double pMomentum = std::sqrt(pEnergy*pEnergy - eMass*eMass);
841 eEnergy = kinEnergy - particleFinalEnergy - pEnergy;
842 G4double eMomentum = std::sqrt(eEnergy*eEnergy - eMass*eMass);
843
845 G4double B3 = -2.*particleMomentum*particleFinalMomentum;
846 G4double cost3interval =
G4Log((A3 + B3*Cmax)/(A3 + B3*Cmin))/B3;
848 G4double cost3 = A3*(expanCost3r6 - 1.)/B3 + Cmin*expanCost3r6;
849 G4double sint3 = std::sqrt((1. - cost3)*(1. + cost3));
850
851 partDirection.
set(sint3*cosp3, sint3*sinp3, cost3);
852
854 muFinalMomentumVector.
set(particleFinalMomentum*sint3*cosp3,
855 particleFinalMomentum*sint3*sinp3,
856 particleFinalMomentum*cost3);
857
858 G4LorentzVector muFourMomentum(particleMomentum, particleMomentumVector);
859 G4LorentzVector muFinalFourMomentum(particleFinalEnergy, muFinalMomentumVector);
862 2.*pEnergy*(kinEnergy - particleFinalEnergy) + 2.*particleMomentumVector[2]*pMomentum -
863 2.*particleFinalMomentum*pMomentum*cost3;
864 G4double B6 = -2.*particleFinalMomentum*pMomentum*(sint3*cosp3*cosp6 + sint3*sinp3*sinp6);
869 G4double t6interval =
G4Log((absA6 + absB6*maxt6)/(absA6 + absB6*mint6))/absB6;
874
875 pDirection.
set(sint6*cosp6, sint6*sinp6, cost6);
877
878 G4ThreeVector auxVec2 = particleMomentumVector - muFinalMomentumVector - pMomentumVector;
880 G4double Bp = particleFinalMomentum*(sint3*cosp3*cosp5 + sint3*sinp3*sinp5) +
881 pMomentum*(sint6*cosp6*cosp5 + sint6*sinp6*sinp5);
882 G4double Cp = -particleMomentum + particleFinalMomentum*cost3 + pMomentum*cost6;
883 G4double A5 = p1mp3mp62 + eMomentum*eMomentum;
890 G4double t5interval = (1./(absA5C5 + absB5*mint5) - 1./(absA5C5 + absB5*maxt5))/absB5;
891 G4double t5 = (-absA5C5 + 1./(1./(absA5C5 + absB5*mint5) - absB5*t5interval*
randNumbs[8]))/absB5;
894
895 eDirection.set(sint5*cosp5, sint5*sinp5, cost5);
896 }
897
898
899
900
901
902
903
904
905
906 auto aParticle1 = new G4DynamicParticle(theElectron, eDirection, eEnergy);
907 auto aParticle2 = new G4DynamicParticle(thePositron, pDirection, pEnergy);
908
909
910 vdp->push_back(aParticle1);
911 vdp->push_back(aParticle2);
912
913
914
918 auto newdp =
new G4DynamicParticle(
particle, muF);
919 vdp->push_back(newdp);
920 } else {
924 }
925
926}
G4double C(G4double temp)
G4double B(G4double temperature)
CLHEP::HepLorentzVector G4LorentzVector
CLHEP::Hep3Vector G4ThreeVector
double dot(const Hep3Vector &) const
void set(double x, double y, double z)
HepLorentzVector & boost(double, double, double)
virtual void flatArray(const int size, double *vect)=0
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
G4double GetTotalEnergy() const
G4ThreeVector GetMomentum() const
G4double GetTotalMomentum() const
G4double FindScaledEnergy(G4int Z, G4double rand, G4double logTkin, G4double yymin, G4double yymax)
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4double SecondaryThreshold() const