| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // GEANT4 Class file
28 //
29 //
30 // File name: G4PolarizedOrePowellAtRestMo
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........oo
51
52 G4PolarizedOrePowellAtRestModel::G4PolarizedOr
53 : G4VPositronAtRestModel("OrePowellPolarized
54 {
55 // DEBUG
56 //G4cout << "G4PolarizedOrePowellAtRestModel
57 }
58
59 //....oooOO0OOooo........oooOO0OOooo........oo
60
61 void G4PolarizedOrePowellAtRestModel::SampleSe
62 std::vector<G4DynamicParticle*>& secParticle
63 {
64 const G4double PositronMass = CLHEP::electro
65 const G4double ymax = 22.0;
66 const G4double sq2 = std::sqrt(2.);
67
68 CLHEP::HepRandomEngine* rndmEngine = G4Rando
69
70 // OrthoPositronium polarization
71
72 // Random mState in 1:1:1 proportion, must b
73 G4int mPos = static_cast<G4int>(std::floor(
74
75 const G4complex iComplex(0., 1.);
76
77 // The u vector defined in global coordinate
78 //quantization along z - axis, must be a use
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, Photon
105 G4ThreeVector PhotonPol1, PhotonPol2, Photon
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 norma
116 r1 = rndmv2[0];
117 r2 = rndmv2[1];
118
119 // energy conservation, with positronium
120 r3 = 2.0 - (r1 + r2);
121 // cosine of angles between photons, fro
122 cos12 = (r3 * r3 - r1 * r1 - r2 * r2) /
123 cos13 = (r2 * r2 - r1 * r1 - r3 * r3) /
124 // request both cosines < 1.
125 } while (std::abs(cos12) > 1 || std::abs(c
126
127 G4double sin12 = std::sqrt((1 + cos12)*(1
128 G4double sin13 = -std::sqrt((1 + cos13)*(1
129
130 // Photon directions in the decay plane (y
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 d
138 G4ThreeVector PhotonPerp1(0., 1., 0.);
139 G4ThreeVector PhotonPerp2(0., cos12, -sin1
140 G4ThreeVector PhotonPerp3(0., cos13, -sin1
141
142 G4ThreeVector xDir(1., 0., 0.);
143
144 // photon polarization vector (a_i in Ore
145 G4double eta1 = CLHEP::pi * G4UniformRand(
146 PhotonPol1 = std::cos(eta1) * PhotonPerp1
147
148 G4double eta2 = CLHEP::pi * G4UniformRand(
149 PhotonPol2 = std::cos(eta2) * PhotonPerp2
150
151 G4double eta3 = CLHEP::pi * G4UniformRand(
152 PhotonPol3 = std::cos(eta3) * PhotonPerp3
153
154 // direction perpendicular to the photon d
155 G4ThreeVector PhotonPrime1 = PhotonPol1.cr
156 G4ThreeVector PhotonPrime2 = PhotonPol2.cr
157 G4ThreeVector PhotonPrime3 = PhotonPol3.cr
158
159 // transition matrix elements (t_i in Ore
160 G4ThreeVector t1(-PhotonPol3 * (PhotonPol1
161 + PhotonPol1 * (PhotonPri
162 - PhotonPrime2 * (PhotonP
163 - PhotonPrime3 * (PhotonP
164 G4ThreeVector t2(-PhotonPol1 * (PhotonPol2
165 + PhotonPol2 * (PhotonPri
166 - PhotonPrime3 * (PhotonP
167 - PhotonPrime1 * (PhotonP
168 G4ThreeVector t3(-PhotonPol2 * (PhotonPol3
169 + PhotonPol3 * (PhotonPri
170 - PhotonPrime1 * (PhotonP
171 - PhotonPrime2 * (PhotonP
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 coordinat
186
187 G4complex H = t(0) * ux + t(1) * uy + t(2)
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
201 aGamma1->SetPolarization(PhotonPol1);
202 secParticles.push_back(aGamma1);
203
204 auto aGamma2 = new G4DynamicParticle(G4Gamma
205 aGamma2->SetPolarization(PhotonPol2);
206 secParticles.push_back(aGamma2);
207
208 auto aGamma3 = new G4DynamicParticle(G4Gamma
209 aGamma3->SetPolarization(PhotonPol3);
210 secParticles.push_back(aGamma3);
211 }
212
213 //....oooOO0OOooo........oooOO0OOooo........oo
214
215 void G4PolarizedOrePowellAtRestModel::PrintGen
216 {
217 G4cout << "Polarized Ore Powell AtRest posit
218 }
219
220 //....oooOO0OOooo........oooOO0OOooo........oo
221