| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // GEANT4 Class file
28 //
29 //
30 // File name: G4PolarizedOrePowellAtRestModel
31 //
32 // Author: I.Semeniouk & D.Bernard
33 //
34 // Creation date: 26 July 2024
35 //
36 // -------------------------------------------------------------------
37 //
38
39 #include "G4PolarizedOrePowellAtRestModel.hh"
40
41 #include "G4DynamicParticle.hh"
42 #include "G4Gamma.hh"
43 #include "G4Material.hh"
44 #include "G4PhysicalConstants.hh"
45 #include "G4RandomDirection.hh"
46 #include "G4RotationMatrix.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "G4ThreeVector.hh"
49
50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
51
52 G4PolarizedOrePowellAtRestModel::G4PolarizedOrePowellAtRestModel()
53 : G4VPositronAtRestModel("OrePowellPolarized")
54 {
55 // DEBUG
56 //G4cout << "G4PolarizedOrePowellAtRestModel in contractor" << G4endl;
57 }
58
59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
60
61 void G4PolarizedOrePowellAtRestModel::SampleSecondaries(
62 std::vector<G4DynamicParticle*>& secParticles, G4double&, const G4Material*) const
63 {
64 const G4double PositronMass = CLHEP::electron_mass_c2;
65 const G4double ymax = 22.0;
66 const G4double sq2 = std::sqrt(2.);
67
68 CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
69
70 // OrthoPositronium polarization
71
72 // Random mState in 1:1:1 proportion, must be a user parameter
73 G4int mPos = static_cast<G4int>(std::floor(3.0 * rndmEngine->flat()) - 1.0);
74
75 const G4complex iComplex(0., 1.);
76
77 // The u vector defined in global coordinate system
78 //quantization along z - axis, must be a user parameter
79 G4complex ux, uy, uz;
80
81 if (mPos == 0) {
82 uz = 1.;
83 uy = 0.;
84 ux = 0.;
85 }
86 else if (mPos == 1) {
87 uz = 0.;
88 uy = iComplex / sq2;
89 ux = 1. / sq2;
90 }
91 else if (mPos == -1) {
92 uz = 0.;
93 uy = -iComplex / sq2;
94 ux = 1. / sq2;
95 }
96
97 G4double r1, r2, r3;
98 G4double cos12, cos13;
99 G4double pdf;
100
101 G4double rndmv2[2];
102 G4double rndmv1;
103
104 G4ThreeVector PhotonDir1, PhotonDir2, PhotonDir3;
105 G4ThreeVector PhotonPol1, PhotonPol2, PhotonPol3;
106
107 // rotate from local to world coordinate;
108 G4RotationMatrix LtoW;
109
110 do {
111 rndmv1 = rndmEngine->flat();
112
113 do {
114 rndmEngine->flatArray(2, rndmv2);
115 // energies of photon1 and photon2 normalized to electron rest mass
116 r1 = rndmv2[0];
117 r2 = rndmv2[1];
118
119 // energy conservation, with positronium assumed = 2 * electron rest mass
120 r3 = 2.0 - (r1 + r2);
121 // cosine of angles between photons, from momentum conservation
122 cos12 = (r3 * r3 - r1 * r1 - r2 * r2) / (2 * r1 * r2);
123 cos13 = (r2 * r2 - r1 * r1 - r3 * r3) / (2 * r1 * r3);
124 // request both cosines < 1.
125 } while (std::abs(cos12) > 1 || std::abs(cos13) > 1);
126
127 G4double sin12 = std::sqrt((1 + cos12)*(1 - cos12));
128 G4double sin13 = -std::sqrt((1 + cos13)*(1 - cos13));
129
130 // Photon directions in the decay plane (y-z), photon 1 along z
131 PhotonDir1 = {0., 0., 1.};
132 PhotonDir2 = {0., sin12, cos12};
133 PhotonDir3 = {0., sin13, cos13};
134
135 // Polarization
136
137 // direction perpendicular to the photon direction, still in the decay plane
138 G4ThreeVector PhotonPerp1(0., 1., 0.);
139 G4ThreeVector PhotonPerp2(0., cos12, -sin12);
140 G4ThreeVector PhotonPerp3(0., cos13, -sin13);
141
142 G4ThreeVector xDir(1., 0., 0.);
143
144 // photon polarization vector (a_i in Ore & Powell)
145 G4double eta1 = CLHEP::pi * G4UniformRand();
146 PhotonPol1 = std::cos(eta1) * PhotonPerp1 + std::sin(eta1) * xDir;
147
148 G4double eta2 = CLHEP::pi * G4UniformRand();
149 PhotonPol2 = std::cos(eta2) * PhotonPerp2 + std::sin(eta2) * xDir;
150
151 G4double eta3 = CLHEP::pi * G4UniformRand();
152 PhotonPol3 = std::cos(eta3) * PhotonPerp3 + std::sin(eta3) * xDir;
153
154 // direction perpendicular to the photon direction and to polarization (a'_i in Ore & Powell)
155 G4ThreeVector PhotonPrime1 = PhotonPol1.cross(PhotonDir1);
156 G4ThreeVector PhotonPrime2 = PhotonPol2.cross(PhotonDir2);
157 G4ThreeVector PhotonPrime3 = PhotonPol3.cross(PhotonDir3);
158
159 // transition matrix elements (t_i in Ore & Powell)
160 G4ThreeVector t1(-PhotonPol3 * (PhotonPol1.dot(PhotonPol2))
161 + PhotonPol1 * (PhotonPrime2.dot(PhotonPrime3))
162 - PhotonPrime2 * (PhotonPrime3.dot(PhotonPol1))
163 - PhotonPrime3 * (PhotonPol1.dot(PhotonPrime2)));
164 G4ThreeVector t2(-PhotonPol1 * (PhotonPol2.dot(PhotonPol3))
165 + PhotonPol2 * (PhotonPrime3.dot(PhotonPrime1))
166 - PhotonPrime3 * (PhotonPrime1.dot(PhotonPol2))
167 - PhotonPrime1 * (PhotonPol2.dot(PhotonPrime3)));
168 G4ThreeVector t3(-PhotonPol2 * (PhotonPol3.dot(PhotonPol1))
169 + PhotonPol3 * (PhotonPrime1.dot(PhotonPrime2))
170 - PhotonPrime1 * (PhotonPrime2.dot(PhotonPol3))
171 - PhotonPrime2 * (PhotonPol3.dot(PhotonPrime1)));
172
173 // Transformation matrix
174
175 G4ThreeVector z = G4RandomDirection();
176 auto tmp = G4RandomDirection();
177 G4ThreeVector x = z.cross(tmp).unit();
178 G4ThreeVector y = z.cross(x);
179
180 LtoW = G4RotationMatrix(x, y, z);
181
182 // G4cout << x << y << z << G4endl;
183
184 G4ThreeVector t = t1 + t2 + t3;
185 t.transform(LtoW); // t in World coordinate system
186
187 G4complex H = t(0) * ux + t(1) * uy + t(2) * uz;
188 pdf = std::abs(conj(H) * H);
189 } while (pdf < ymax * rndmv1);
190 // END of Sampling
191
192 PhotonDir1.transform(LtoW);
193 PhotonDir2.transform(LtoW);
194 PhotonDir3.transform(LtoW);
195
196 PhotonPol1.transform(LtoW);
197 PhotonPol2.transform(LtoW);
198 PhotonPol3.transform(LtoW);
199
200 auto aGamma1 = new G4DynamicParticle(G4Gamma::Gamma(), PhotonDir1, r1 * PositronMass);
201 aGamma1->SetPolarization(PhotonPol1);
202 secParticles.push_back(aGamma1);
203
204 auto aGamma2 = new G4DynamicParticle(G4Gamma::Gamma(), PhotonDir2, r2 * PositronMass);
205 aGamma2->SetPolarization(PhotonPol2);
206 secParticles.push_back(aGamma2);
207
208 auto aGamma3 = new G4DynamicParticle(G4Gamma::Gamma(), PhotonDir3, r3 * PositronMass);
209 aGamma3->SetPolarization(PhotonPol3);
210 secParticles.push_back(aGamma3);
211 }
212
213 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
214
215 void G4PolarizedOrePowellAtRestModel::PrintGeneratorInformation() const
216 {
217 G4cout << "Polarized Ore Powell AtRest positron 3-gamma annihilation model" << G4endl;
218 }
219
220 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
221