58{
59 static const G4double PositronMass = CLHEP::electron_mass_c2;
61
62 CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
63
72
75
76 do {
77 rndmv1 = rndmEngine->
flat();
78
79 do {
81
82 r1 = rndmv2[0];
83 r2 = rndmv2[1];
84
85 r3 = 2.0 - (r1+r2);
86
87 cos12=(r3*r3 - r1*r1 -r2*r2)/(2*r1*r2);
88 cos13=(r2*r2 - r1*r1 -r3*r3)/(2*r1*r3);
89
90 } while ( std::abs(cos12) > 1 || std::abs(cos13) > 1 );
91
92 sin12 = std::sqrt((1 + cos12)*(1 - cos12));
93 sin13 = -std::sqrt((1 + cos13)*(1 - cos13));
94
95 G4double cos23=cos12*cos13+sin12*sin13;
96
97 pdf = (1 - cos12)*(1 - cos12) + (1 - cos13)*(1 - cos13) + (1 - cos23)*(1 - cos23);
98 } while ( pdf < ymax * rndmv1 );
99
100
101
105
106
108
109 PhotonMomentum1.rotateUz(dir1);
110 PhotonMomentum2.rotateUz(dir1);
111 PhotonMomentum3.rotateUz(dir1);
112
113 auto aGamma1 =
new G4DynamicParticle(
G4Gamma::Gamma(), PhotonMomentum1,
114 r1 * PositronMass);
115
116
119 pol1.rotateUz(PhotonMomentum1);
120 aGamma1->SetPolarization(pol1);
121 secParticles.push_back(aGamma1);
122
123 auto aGamma2 =
new G4DynamicParticle(
G4Gamma::Gamma(), PhotonMomentum2,
124 r2 * PositronMass);
125
128 pol2.rotateUz(PhotonMomentum2);
129 aGamma2->SetPolarization(pol2);
130 secParticles.push_back(aGamma2);
131
132 auto aGamma3 =
new G4DynamicParticle(
G4Gamma::Gamma(), PhotonMomentum3,
133 r3 * PositronMass);
134
137 pol3.rotateUz(PhotonMomentum3);
138 aGamma3->SetPolarization(pol3);
139 secParticles.push_back(aGamma3);
140
141}
G4ThreeVector G4RandomDirection()
CLHEP::Hep3Vector G4ThreeVector
virtual void flatArray(const int size, double *vect)=0