| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // ------------------------------------------- 26 // -------------------------------------------------------------------
27 // 27 //
28 // Geant4 Class file 28 // Geant4 Class file
29 // 29 //
30 // File name: G4PolarizedComptonModel 30 // File name: G4PolarizedComptonModel
31 // 31 //
32 // Author: Andreas Schaelicke 32 // Author: Andreas Schaelicke
33 33
34 #include "G4PolarizedComptonModel.hh" 34 #include "G4PolarizedComptonModel.hh"
35 35
36 #include "G4Exp.hh" 36 #include "G4Exp.hh"
37 #include "G4Log.hh" 37 #include "G4Log.hh"
38 #include "G4ParticleChangeForGamma.hh" 38 #include "G4ParticleChangeForGamma.hh"
39 #include "G4PhysicalConstants.hh" 39 #include "G4PhysicalConstants.hh"
40 #include "G4PolarizationManager.hh" 40 #include "G4PolarizationManager.hh"
41 #include "G4PolarizationHelper.hh" 41 #include "G4PolarizationHelper.hh"
42 #include "G4PolarizedComptonXS.hh" 42 #include "G4PolarizedComptonXS.hh"
43 #include "G4StokesVector.hh" 43 #include "G4StokesVector.hh"
44 #include "G4SystemOfUnits.hh" 44 #include "G4SystemOfUnits.hh"
45 45
46 //....oooOO0OOooo........oooOO0OOooo........oo 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47 G4PolarizedComptonModel::G4PolarizedComptonMod 47 G4PolarizedComptonModel::G4PolarizedComptonModel(const G4ParticleDefinition*,
48 48 const G4String& nam)
49 : G4KleinNishinaCompton(nullptr, nam) 49 : G4KleinNishinaCompton(nullptr, nam)
50 , fVerboseLevel(0) 50 , fVerboseLevel(0)
51 { 51 {
52 fCrossSectionCalculator = new G4PolarizedCom 52 fCrossSectionCalculator = new G4PolarizedComptonXS();
53 fBeamPolarization = G4StokesVector::ZE 53 fBeamPolarization = G4StokesVector::ZERO;
54 fTargetPolarization = G4StokesVector::ZE 54 fTargetPolarization = G4StokesVector::ZERO;
55 } 55 }
56 56
57 //....oooOO0OOooo........oooOO0OOooo........oo 57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
58 G4PolarizedComptonModel::~G4PolarizedComptonMo 58 G4PolarizedComptonModel::~G4PolarizedComptonModel()
59 { 59 {
60 delete fCrossSectionCalculator; 60 delete fCrossSectionCalculator;
61 } 61 }
62 62
63 //....oooOO0OOooo........oooOO0OOooo........oo 63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
64 G4double G4PolarizedComptonModel::ComputeAsymm 64 G4double G4PolarizedComptonModel::ComputeAsymmetryPerAtom(G4double gammaEnergy,
65 65 G4double /*Z*/)
66 { 66 {
67 G4double asymmetry = 0.0; 67 G4double asymmetry = 0.0;
68 68
69 G4double k0 = gammaEnergy / electron_mass_c2 69 G4double k0 = gammaEnergy / electron_mass_c2;
70 G4double k1 = 1. + 2. * k0; 70 G4double k1 = 1. + 2. * k0;
71 71
72 asymmetry = -k0; 72 asymmetry = -k0;
73 asymmetry *= 73 asymmetry *=
74 (k0 + 1.) * sqr(k1) * G4Log(k1) - 2. * k0 74 (k0 + 1.) * sqr(k1) * G4Log(k1) - 2. * k0 * (5. * sqr(k0) + 4. * k0 + 1.);
75 asymmetry /= ((k0 - 2.) * k0 - 2.) * sqr(k1) 75 asymmetry /= ((k0 - 2.) * k0 - 2.) * sqr(k1) * G4Log(k1) +
76 2. * k0 * (k0 * (k0 + 1.) * (k0 76 2. * k0 * (k0 * (k0 + 1.) * (k0 + 8.) + 2.);
77 77
78 if(asymmetry > 1.) 78 if(asymmetry > 1.)
79 { 79 {
80 G4ExceptionDescription ed; 80 G4ExceptionDescription ed;
81 ed << "ERROR in G4PolarizedComptonModel::C 81 ed << "ERROR in G4PolarizedComptonModel::ComputeAsymmetryPerAtom.\n"
82 << " asymmetry = " << asymmetry << "\n" 82 << " asymmetry = " << asymmetry << "\n";
83 G4Exception("G4PolarizedComptonModel::Comp 83 G4Exception("G4PolarizedComptonModel::ComputeAsymmetryPerAtom", "pol035",
84 JustWarning, ed); 84 JustWarning, ed);
85 } 85 }
86 86
87 return asymmetry; 87 return asymmetry;
88 } 88 }
89 89
90 //....oooOO0OOooo........oooOO0OOooo........oo 90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
91 G4double G4PolarizedComptonModel::ComputeCross 91 G4double G4PolarizedComptonModel::ComputeCrossSectionPerAtom(
92 const G4ParticleDefinition* pd, G4double kin 92 const G4ParticleDefinition* pd, G4double kinEnergy, G4double Z, G4double A,
93 G4double cut, G4double emax) 93 G4double cut, G4double emax)
94 { 94 {
95 G4double xs = G4KleinNishinaCompton::Compute 95 G4double xs = G4KleinNishinaCompton::ComputeCrossSectionPerAtom(
96 pd, kinEnergy, Z, A, cut, emax); 96 pd, kinEnergy, Z, A, cut, emax);
97 G4double polzz = fBeamPolarization.p3() * fT 97 G4double polzz = fBeamPolarization.p3() * fTargetPolarization.z();
98 if(polzz > 0.0) 98 if(polzz > 0.0)
99 { 99 {
100 G4double asym = ComputeAsymmetryPerAtom(ki 100 G4double asym = ComputeAsymmetryPerAtom(kinEnergy, Z);
101 xs *= (1. + polzz * asym); 101 xs *= (1. + polzz * asym);
102 } 102 }
103 return xs; 103 return xs;
104 } 104 }
105 105
106 //....oooOO0OOooo........oooOO0OOooo........oo 106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
107 void G4PolarizedComptonModel::SampleSecondarie 107 void G4PolarizedComptonModel::SampleSecondaries(
108 std::vector<G4DynamicParticle*>* fvect, cons 108 std::vector<G4DynamicParticle*>* fvect, const G4MaterialCutsCouple*,
109 const G4DynamicParticle* aDynamicGamma, G4do 109 const G4DynamicParticle* aDynamicGamma, G4double, G4double)
110 { 110 {
111 // do nothing below the threshold 111 // do nothing below the threshold
112 if(aDynamicGamma->GetKineticEnergy() <= LowE 112 if(aDynamicGamma->GetKineticEnergy() <= LowEnergyLimit())
113 { 113 {
114 return; 114 return;
115 } 115 }
116 116
117 const G4Track* aTrack = fParticleChang 117 const G4Track* aTrack = fParticleChange->GetCurrentTrack();
118 G4VPhysicalVolume* aPVolume = aTrack->GetVol 118 G4VPhysicalVolume* aPVolume = aTrack->GetVolume();
119 G4LogicalVolume* aLVolume = aPVolume->GetL 119 G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
120 120
121 if(fVerboseLevel >= 1) 121 if(fVerboseLevel >= 1)
122 { 122 {
123 G4cout << "G4PolarizedComptonModel::Sample 123 G4cout << "G4PolarizedComptonModel::SampleSecondaries in "
124 << aLVolume->GetName() << G4endl; 124 << aLVolume->GetName() << G4endl;
125 } 125 }
126 G4PolarizationManager* polarizationManager = 126 G4PolarizationManager* polarizationManager =
127 G4PolarizationManager::GetInstance(); 127 G4PolarizationManager::GetInstance();
128 128
129 // obtain polarization of the beam 129 // obtain polarization of the beam
130 fBeamPolarization = G4StokesVector(aDynamicG 130 fBeamPolarization = G4StokesVector(aDynamicGamma->GetPolarization());
131 fBeamPolarization.SetPhoton(); 131 fBeamPolarization.SetPhoton();
132 132
133 // obtain polarization of the media 133 // obtain polarization of the media
134 G4bool targetIsPolarized = polarizationManag 134 G4bool targetIsPolarized = polarizationManager->IsPolarized(aLVolume);
135 fTargetPolarization = polarizationManager->G 135 fTargetPolarization = polarizationManager->GetVolumePolarization(aLVolume);
136 136
137 // if beam is linear polarized or target is 137 // if beam is linear polarized or target is transversely polarized
138 // determine the angle to x-axis 138 // determine the angle to x-axis
139 // (assumes same PRF as in the polarization 139 // (assumes same PRF as in the polarization definition)
140 G4ThreeVector gamDirection0 = aDynamicGamma- 140 G4ThreeVector gamDirection0 = aDynamicGamma->GetMomentumDirection();
141 141
142 // transfer fTargetPolarization 142 // transfer fTargetPolarization
143 // into the gamma frame (problem electron is 143 // into the gamma frame (problem electron is at rest)
144 if(targetIsPolarized) 144 if(targetIsPolarized)
145 { 145 {
146 fTargetPolarization.rotateUz(gamDirection0 146 fTargetPolarization.rotateUz(gamDirection0);
147 } 147 }
148 // The scattered gamma energy is sampled acc 148 // The scattered gamma energy is sampled according to
149 // Klein - Nishina formula. 149 // Klein - Nishina formula.
150 // The random number techniques of Butcher & 150 // The random number techniques of Butcher & Messel are used
151 // (Nuc Phys 20(1960),15). 151 // (Nuc Phys 20(1960),15).
152 // Note : Effects due to binding of atomic e 152 // Note : Effects due to binding of atomic electrons are neglected.
153 153
154 G4double gamEnergy0 = aDynamicGamma->GetKine 154 G4double gamEnergy0 = aDynamicGamma->GetKineticEnergy();
155 G4double E0_m = gamEnergy0 / electron_ 155 G4double E0_m = gamEnergy0 / electron_mass_c2;
156 156
157 // sample the energy rate of the scattered g 157 // sample the energy rate of the scattered gamma
158 G4double epsilon, sint2; 158 G4double epsilon, sint2;
159 G4double onecost = 0.0; 159 G4double onecost = 0.0;
160 G4double Phi = 0.0; 160 G4double Phi = 0.0;
161 G4double greject = 1.0; 161 G4double greject = 1.0;
162 G4double cosTeta = 1.0; 162 G4double cosTeta = 1.0;
163 G4double sinTeta = 0.0; 163 G4double sinTeta = 0.0;
164 164
165 G4double eps0 = 1. / (1. + 2. * E0_m); 165 G4double eps0 = 1. / (1. + 2. * E0_m);
166 G4double epsilon0sq = eps0 * eps0; 166 G4double epsilon0sq = eps0 * eps0;
167 G4double alpha1 = -G4Log(eps0); 167 G4double alpha1 = -G4Log(eps0);
168 G4double alpha2 = alpha1 + 0.5 * (1. - e 168 G4double alpha2 = alpha1 + 0.5 * (1. - epsilon0sq);
169 169
170 G4double polarization = fBeamPolarization.p3 170 G4double polarization = fBeamPolarization.p3() * fTargetPolarization.p3();
171 171
172 CLHEP::HepRandomEngine* rndmEngineMod = G4Ra 172 CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine();
173 G4int nloop = 0; 173 G4int nloop = 0;
174 G4bool end = fals 174 G4bool end = false;
175 175
176 G4double rndm[3]; 176 G4double rndm[3];
177 177
178 do 178 do
179 { 179 {
180 do 180 do
181 { 181 {
182 ++nloop; 182 ++nloop;
183 // false interaction if too many iterati 183 // false interaction if too many iterations
184 if(nloop > fLoopLim) 184 if(nloop > fLoopLim)
185 { 185 {
186 PrintWarning(aDynamicGamma, nloop, gre 186 PrintWarning(aDynamicGamma, nloop, greject, onecost, Phi,
187 "too many iterations"); 187 "too many iterations");
188 return; 188 return;
189 } 189 }
190 190
191 // 3 random numbers to sample scattering 191 // 3 random numbers to sample scattering
192 rndmEngineMod->flatArray(3, rndm); 192 rndmEngineMod->flatArray(3, rndm);
193 193
194 if(alpha1 > alpha2 * rndm[0]) 194 if(alpha1 > alpha2 * rndm[0])
195 { 195 {
196 epsilon = G4Exp(-alpha1 * rndm[1]); 196 epsilon = G4Exp(-alpha1 * rndm[1]);
197 } 197 }
198 else 198 else
199 { 199 {
200 epsilon = std::sqrt(epsilon0sq + (1. - 200 epsilon = std::sqrt(epsilon0sq + (1. - epsilon0sq) * rndm[1]);
201 } 201 }
202 202
203 onecost = (1. - epsilon) / (epsilon * E0 203 onecost = (1. - epsilon) / (epsilon * E0_m);
204 sint2 = onecost * (2. - onecost); 204 sint2 = onecost * (2. - onecost);
205 205
206 G4double gdiced = 2. * (1. / epsilon + e 206 G4double gdiced = 2. * (1. / epsilon + epsilon);
207 G4double gdist = 1. / epsilon + epsilon 207 G4double gdist = 1. / epsilon + epsilon - sint2 -
208 polarization * (1. / ep 208 polarization * (1. / epsilon - epsilon) * (1. - onecost);
209 209
210 greject = gdist / gdiced; 210 greject = gdist / gdiced;
211 211
212 if(greject > 1.0) 212 if(greject > 1.0)
213 { 213 {
214 PrintWarning(aDynamicGamma, nloop, gre 214 PrintWarning(aDynamicGamma, nloop, greject, onecost, Phi,
215 "theta majoranta wrong"); 215 "theta majoranta wrong");
216 } 216 }
217 // Loop checking, 03-Aug-2015, Vladimir 217 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
218 } while(greject < rndm[2]); 218 } while(greject < rndm[2]);
219 219
220 // assuming phi loop successful 220 // assuming phi loop successful
221 end = true; 221 end = true;
222 222
223 // scattered gamma angles. ( Z - axis alon 223 // scattered gamma angles. ( Z - axis along the parent gamma)
224 cosTeta = 1. - onecost; 224 cosTeta = 1. - onecost;
225 sinTeta = std::sqrt(sint2); 225 sinTeta = std::sqrt(sint2);
226 do 226 do
227 { 227 {
228 ++nloop; 228 ++nloop;
229 229
230 // 2 random numbers to sample scattering 230 // 2 random numbers to sample scattering
231 rndmEngineMod->flatArray(2, rndm); 231 rndmEngineMod->flatArray(2, rndm);
232 232
233 // false interaction if too many iterati 233 // false interaction if too many iterations
234 Phi = twopi * rndm[0]; 234 Phi = twopi * rndm[0];
235 if(nloop > fLoopLim) 235 if(nloop > fLoopLim)
236 { 236 {
237 PrintWarning(aDynamicGamma, nloop, gre 237 PrintWarning(aDynamicGamma, nloop, greject, onecost, Phi,
238 "too many iterations"); 238 "too many iterations");
239 return; 239 return;
240 } 240 }
241 241
242 G4double gdiced = 1. / epsilon + epsilon 242 G4double gdiced = 1. / epsilon + epsilon - sint2 +
243 std::abs(fBeamPolariza 243 std::abs(fBeamPolarization.p3()) *
244 (std::abs((1. / epsi 244 (std::abs((1. / epsilon - epsilon) * cosTeta *
245 fTargetPol 245 fTargetPolarization.p3()) +
246 (1. - epsilon) * si 246 (1. - epsilon) * sinTeta *
247 (std::sqrt(sqr(fT 247 (std::sqrt(sqr(fTargetPolarization.p1()) +
248 sqr(fT 248 sqr(fTargetPolarization.p2())))) +
249 sint2 * (std::sqrt(sqr 249 sint2 * (std::sqrt(sqr(fBeamPolarization.p1()) +
250 sqr 250 sqr(fBeamPolarization.p2())));
251 251
252 G4double gdist = 252 G4double gdist =
253 1. / epsilon + epsilon - sint2 + 253 1. / epsilon + epsilon - sint2 +
254 fBeamPolarization.p3() * 254 fBeamPolarization.p3() *
255 ((1. / epsilon - epsilon) * cosTeta 255 ((1. / epsilon - epsilon) * cosTeta * fTargetPolarization.p3() +
256 (1. - epsilon) * sinTeta * 256 (1. - epsilon) * sinTeta *
257 (std::cos(Phi) * fTargetPolarizat 257 (std::cos(Phi) * fTargetPolarization.p1() +
258 std::sin(Phi) * fTargetPolarizat 258 std::sin(Phi) * fTargetPolarization.p2())) -
259 sint2 * (std::cos(2. * Phi) * fBeamPol 259 sint2 * (std::cos(2. * Phi) * fBeamPolarization.p1() +
260 std::sin(2. * Phi) * fBeamPol 260 std::sin(2. * Phi) * fBeamPolarization.p2());
261 greject = gdist / gdiced; 261 greject = gdist / gdiced;
262 262
263 if(greject > 1.0) 263 if(greject > 1.0)
264 { 264 {
265 PrintWarning(aDynamicGamma, nloop, gre 265 PrintWarning(aDynamicGamma, nloop, greject, onecost, Phi,
266 "phi majoranta wrong"); 266 "phi majoranta wrong");
267 } 267 }
268 268
269 if(greject < 1.e-3) 269 if(greject < 1.e-3)
270 { 270 {
271 PrintWarning(aDynamicGamma, nloop, gre 271 PrintWarning(aDynamicGamma, nloop, greject, onecost, Phi,
272 "phi loop ineffective"); 272 "phi loop ineffective");
273 // restart theta loop 273 // restart theta loop
274 end = false; 274 end = false;
275 break; 275 break;
276 } 276 }
277 277
278 // Loop checking, 03-Aug-2015, Vladimir 278 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
279 } while(greject < rndm[1]); 279 } while(greject < rndm[1]);
280 } while(!end); 280 } while(!end);
281 G4double dirx = sinTeta * std::cos(Phi); 281 G4double dirx = sinTeta * std::cos(Phi);
282 G4double diry = sinTeta * std::sin(Phi); 282 G4double diry = sinTeta * std::sin(Phi);
283 G4double dirz = cosTeta; 283 G4double dirz = cosTeta;
284 284
285 // update G4VParticleChange for the scattere 285 // update G4VParticleChange for the scattered gamma
286 G4ThreeVector gamDirection1(dirx, diry, dirz 286 G4ThreeVector gamDirection1(dirx, diry, dirz);
287 gamDirection1.rotateUz(gamDirection0); 287 gamDirection1.rotateUz(gamDirection0);
288 G4double gamEnergy1 = epsilon * gamEnergy0; 288 G4double gamEnergy1 = epsilon * gamEnergy0;
289 289
290 G4double edep = 0.0; 290 G4double edep = 0.0;
291 if(gamEnergy1 > lowestSecondaryEnergy) 291 if(gamEnergy1 > lowestSecondaryEnergy)
292 { 292 {
293 fParticleChange->ProposeMomentumDirection( 293 fParticleChange->ProposeMomentumDirection(gamDirection1);
294 fParticleChange->SetProposedKineticEnergy( 294 fParticleChange->SetProposedKineticEnergy(gamEnergy1);
295 } 295 }
296 else 296 else
297 { 297 {
298 fParticleChange->ProposeTrackStatus(fStopA 298 fParticleChange->ProposeTrackStatus(fStopAndKill);
299 fParticleChange->SetProposedKineticEnergy( 299 fParticleChange->SetProposedKineticEnergy(0.0);
300 edep = gamEnergy1; 300 edep = gamEnergy1;
301 } 301 }
302 302
303 // calculate Stokes vector of final state ph 303 // calculate Stokes vector of final state photon and electron
304 G4ThreeVector nInteractionFrame = 304 G4ThreeVector nInteractionFrame =
305 G4PolarizationHelper::GetFrame(gamDirectio 305 G4PolarizationHelper::GetFrame(gamDirection1, gamDirection0);
306 306
307 // transfer fBeamPolarization and fTargetPol 307 // transfer fBeamPolarization and fTargetPolarization
308 // into the interaction frame (note electron 308 // into the interaction frame (note electron is in gamma frame)
309 if(fVerboseLevel >= 1) 309 if(fVerboseLevel >= 1)
310 { 310 {
311 G4cout << "=============================== 311 G4cout << "========================================" << G4endl;
312 G4cout << " nInteractionFrame = " << nInte 312 G4cout << " nInteractionFrame = " << nInteractionFrame << G4endl;
313 G4cout << " GammaDirection0 = " << gamDire 313 G4cout << " GammaDirection0 = " << gamDirection0 << G4endl;
314 G4cout << " gammaPolarization = " << fBeam 314 G4cout << " gammaPolarization = " << fBeamPolarization << G4endl;
315 G4cout << " electronPolarization = " << fT 315 G4cout << " electronPolarization = " << fTargetPolarization << G4endl;
316 } 316 }
317 317
318 fBeamPolarization.InvRotateAz(nInteractionFr 318 fBeamPolarization.InvRotateAz(nInteractionFrame, gamDirection0);
319 fTargetPolarization.InvRotateAz(nInteraction 319 fTargetPolarization.InvRotateAz(nInteractionFrame, gamDirection0);
320 320
321 if(fVerboseLevel >= 1) 321 if(fVerboseLevel >= 1)
322 { 322 {
323 G4cout << "------------------------------- 323 G4cout << "----------------------------------------" << G4endl;
324 G4cout << " gammaPolarization = " << fBeam 324 G4cout << " gammaPolarization = " << fBeamPolarization << G4endl;
325 G4cout << " electronPolarization = " << fT 325 G4cout << " electronPolarization = " << fTargetPolarization << G4endl;
326 G4cout << "------------------------------- 326 G4cout << "----------------------------------------" << G4endl;
327 } 327 }
328 328
329 // initialize the polarization transfer matr 329 // initialize the polarization transfer matrix
330 fCrossSectionCalculator->Initialize(epsilon, 330 fCrossSectionCalculator->Initialize(epsilon, E0_m, 0., fBeamPolarization,
331 fTargetP 331 fTargetPolarization, 2);
332 332
333 if(gamEnergy1 > lowestSecondaryEnergy) 333 if(gamEnergy1 > lowestSecondaryEnergy)
334 { 334 {
335 // in interaction frame 335 // in interaction frame
336 // calculate polarization transfer to the 336 // calculate polarization transfer to the photon (in interaction plane)
337 fFinalGammaPolarization = fCrossSectionCal 337 fFinalGammaPolarization = fCrossSectionCalculator->GetPol2();
338 if(fVerboseLevel >= 1) 338 if(fVerboseLevel >= 1)
339 { 339 {
340 G4cout << " gammaPolarization1 = " << fF 340 G4cout << " gammaPolarization1 = " << fFinalGammaPolarization << G4endl;
341 } 341 }
342 fFinalGammaPolarization.SetPhoton(); 342 fFinalGammaPolarization.SetPhoton();
343 343
344 // translate polarization into particle re 344 // translate polarization into particle reference frame
345 fFinalGammaPolarization.RotateAz(nInteract 345 fFinalGammaPolarization.RotateAz(nInteractionFrame, gamDirection1);
346 if(fFinalGammaPolarization.mag() > 1. + 1. 346 if(fFinalGammaPolarization.mag() > 1. + 1.e-8)
347 { 347 {
348 G4ExceptionDescription ed; 348 G4ExceptionDescription ed;
349 ed << "ERROR in Polarizaed Compton Scatt 349 ed << "ERROR in Polarizaed Compton Scattering !\n";
350 ed << "Polarization of final photon more 350 ed << "Polarization of final photon more than 100%.\n";
351 ed << fFinalGammaPolarization 351 ed << fFinalGammaPolarization
352 << " mag = " << fFinalGammaPolarizati 352 << " mag = " << fFinalGammaPolarization.mag() << "\n";
353 G4Exception("G4PolarizedComptonModel::Sa 353 G4Exception("G4PolarizedComptonModel::SampleSecondaries", "pol033",
354 FatalException, ed); 354 FatalException, ed);
355 } 355 }
356 // store polarization vector 356 // store polarization vector
357 fParticleChange->ProposePolarization(fFina 357 fParticleChange->ProposePolarization(fFinalGammaPolarization);
358 if(fVerboseLevel >= 1) 358 if(fVerboseLevel >= 1)
359 { 359 {
360 G4cout << " gammaPolarization1 = " << fF 360 G4cout << " gammaPolarization1 = " << fFinalGammaPolarization << G4endl;
361 G4cout << " GammaDirection1 = " << gamDi 361 G4cout << " GammaDirection1 = " << gamDirection1 << G4endl;
362 } 362 }
363 } 363 }
364 364
365 // kinematic of the scattered electron 365 // kinematic of the scattered electron
366 G4double eKinEnergy = gamEnergy0 - gamEnergy 366 G4double eKinEnergy = gamEnergy0 - gamEnergy1;
367 367
368 if(eKinEnergy > lowestSecondaryEnergy) 368 if(eKinEnergy > lowestSecondaryEnergy)
369 { 369 {
370 G4ThreeVector eDirection = 370 G4ThreeVector eDirection =
371 gamEnergy0 * gamDirection0 - gamEnergy1 371 gamEnergy0 * gamDirection0 - gamEnergy1 * gamDirection1;
372 eDirection = eDirection.unit(); 372 eDirection = eDirection.unit();
373 373
374 finalElectronPolarization = fCrossSectionC 374 finalElectronPolarization = fCrossSectionCalculator->GetPol3();
375 if(fVerboseLevel >= 1) 375 if(fVerboseLevel >= 1)
376 { 376 {
377 G4cout << " electronPolarization1 = " << 377 G4cout << " electronPolarization1 = " << finalElectronPolarization
378 << G4endl; 378 << G4endl;
379 } 379 }
380 // transfer into particle reference frame 380 // transfer into particle reference frame
381 finalElectronPolarization.RotateAz(nIntera 381 finalElectronPolarization.RotateAz(nInteractionFrame, eDirection);
382 if(fVerboseLevel >= 1) 382 if(fVerboseLevel >= 1)
383 { 383 {
384 G4cout << " electronPolarization1 = " << 384 G4cout << " electronPolarization1 = " << finalElectronPolarization
385 << G4endl << " ElecDirection = " 385 << G4endl << " ElecDirection = " << eDirection << G4endl;
386 } 386 }
387 387
388 // create G4DynamicParticle object for the 388 // create G4DynamicParticle object for the electron.
389 G4DynamicParticle* aElectron = 389 G4DynamicParticle* aElectron =
390 new G4DynamicParticle(theElectron, eDire 390 new G4DynamicParticle(theElectron, eDirection, eKinEnergy);
391 // store polarization vector 391 // store polarization vector
392 if(finalElectronPolarization.mag() > 1. + 392 if(finalElectronPolarization.mag() > 1. + 1.e-8)
393 { 393 {
394 G4ExceptionDescription ed; 394 G4ExceptionDescription ed;
395 ed << "ERROR in Polarized Compton Scatte 395 ed << "ERROR in Polarized Compton Scattering !\n";
396 ed << "Polarization of final electron mo 396 ed << "Polarization of final electron more than 100%.\n";
397 ed << finalElectronPolarization 397 ed << finalElectronPolarization
398 << " mag = " << finalElectronPolariza 398 << " mag = " << finalElectronPolarization.mag() << G4endl;
399 G4Exception("G4PolarizedComptonModel::Sa 399 G4Exception("G4PolarizedComptonModel::SampleSecondaries", "pol034",
400 FatalException, ed); 400 FatalException, ed);
401 } 401 }
402 aElectron->SetPolarization(finalElectronPo 402 aElectron->SetPolarization(finalElectronPolarization.p1(),
403 finalElectronPo 403 finalElectronPolarization.p2(),
404 finalElectronPo 404 finalElectronPolarization.p3());
405 fvect->push_back(aElectron); 405 fvect->push_back(aElectron);
406 } 406 }
407 else 407 else
408 { 408 {
409 edep += eKinEnergy; 409 edep += eKinEnergy;
410 } 410 }
411 // energy balance 411 // energy balance
412 if(edep > 0.0) 412 if(edep > 0.0)
413 { 413 {
414 fParticleChange->ProposeLocalEnergyDeposit 414 fParticleChange->ProposeLocalEnergyDeposit(edep);
415 } 415 }
416 } 416 }
417 417
418 //....oooOO0OOooo........oooOO0OOooo........oo 418 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
419 void G4PolarizedComptonModel::PrintWarning(con 419 void G4PolarizedComptonModel::PrintWarning(const G4DynamicParticle* dp,
420 G4i 420 G4int nloop, G4double grej,
421 G4d 421 G4double onecos, G4double phi,
422 con 422 const G4String sss) const
423 { 423 {
424 G4ExceptionDescription ed; 424 G4ExceptionDescription ed;
425 ed << "Problem of scattering sampling: " << 425 ed << "Problem of scattering sampling: " << sss << "\n"
426 << "Niter= " << nloop << " grej= " << gre 426 << "Niter= " << nloop << " grej= " << grej
427 << " cos(theta)= " << 1.0 - onecos << " p 427 << " cos(theta)= " << 1.0 - onecos << " phi= " << phi << "\n"
428 << "Gamma E(MeV)= " << dp->GetKineticEner 428 << "Gamma E(MeV)= " << dp->GetKineticEnergy() / MeV
429 << " dir= " << dp->GetMomentumDirection() 429 << " dir= " << dp->GetMomentumDirection()
430 << " pol= " << dp->GetPolarization(); 430 << " pol= " << dp->GetPolarization();
431 G4Exception("G4PolarizedComptonModel::Sample 431 G4Exception("G4PolarizedComptonModel::SampleSecondaries", "em0044",
432 JustWarning, ed, ""); 432 JustWarning, ed, "");
433 } 433 }
434 434