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 | // $Id: G4PolarizedComptonModel.cc,v 1.4 2007/05/23 08:52:20 vnivanch Exp $ |
---|
28 | // GEANT4 tag $Name: geant4-09-03 $ |
---|
29 | // |
---|
30 | // ------------------------------------------------------------------- |
---|
31 | // |
---|
32 | // GEANT4 Class file |
---|
33 | // |
---|
34 | // |
---|
35 | // File name: G4PolarizedComptonModel |
---|
36 | // |
---|
37 | // Author: Andreas Schaelicke |
---|
38 | // |
---|
39 | // Creation date: 01.05.2005 |
---|
40 | // |
---|
41 | // Modifications: |
---|
42 | // 18-07-06 use newly calculated cross sections (P. Starovoitov) |
---|
43 | // 21-08-05 update interface (A. Schaelicke) |
---|
44 | // |
---|
45 | // Class Description: |
---|
46 | // |
---|
47 | // ------------------------------------------------------------------- |
---|
48 | // |
---|
49 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
50 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
51 | |
---|
52 | #include "G4PolarizedComptonModel.hh" |
---|
53 | #include "G4Electron.hh" |
---|
54 | #include "G4Gamma.hh" |
---|
55 | #include "Randomize.hh" |
---|
56 | #include "G4DataVector.hh" |
---|
57 | #include "G4ParticleChangeForGamma.hh" |
---|
58 | |
---|
59 | |
---|
60 | #include "G4StokesVector.hh" |
---|
61 | #include "G4PolarizationManager.hh" |
---|
62 | #include "G4PolarizationHelper.hh" |
---|
63 | #include "G4PolarizedComptonCrossSection.hh" |
---|
64 | |
---|
65 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
66 | |
---|
67 | G4PolarizedComptonModel::G4PolarizedComptonModel(const G4ParticleDefinition*, |
---|
68 | const G4String& nam) |
---|
69 | : G4KleinNishinaCompton(0,nam), |
---|
70 | verboseLevel(0) |
---|
71 | { |
---|
72 | crossSectionCalculator=new G4PolarizedComptonCrossSection(); |
---|
73 | } |
---|
74 | |
---|
75 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
76 | |
---|
77 | G4PolarizedComptonModel::~G4PolarizedComptonModel() |
---|
78 | { |
---|
79 | if (crossSectionCalculator) delete crossSectionCalculator; |
---|
80 | } |
---|
81 | |
---|
82 | |
---|
83 | |
---|
84 | G4double G4PolarizedComptonModel::ComputeAsymmetryPerAtom |
---|
85 | (G4double gammaEnergy, G4double /*Z*/) |
---|
86 | |
---|
87 | { |
---|
88 | G4double asymmetry = 0.0 ; |
---|
89 | |
---|
90 | G4double k0 = gammaEnergy / electron_mass_c2 ; |
---|
91 | G4double k1 = 1 + 2*k0 ; |
---|
92 | |
---|
93 | asymmetry = -k0; |
---|
94 | asymmetry *= (k0 + 1.)*sqr(k1)*std::log(k1) - 2.*k0*(5.*sqr(k0) + 4.*k0 + 1.); |
---|
95 | asymmetry /= ((k0 - 2.)*k0 -2.)*sqr(k1)*std::log(k1) + 2.*k0*(k0*(k0 + 1.)*(k0 + 8.) + 2.); |
---|
96 | |
---|
97 | // G4cout<<"energy = "<<GammaEnergy<<" asymmetry = "<<asymmetry<<"\t\t GAM = "<<k0<<G4endl; |
---|
98 | if (asymmetry>1.) G4cout<<"ERROR in G4PolarizedComptonModel::ComputeAsymmetryPerAtom"<<G4endl; |
---|
99 | |
---|
100 | return asymmetry; |
---|
101 | } |
---|
102 | |
---|
103 | |
---|
104 | G4double G4PolarizedComptonModel::ComputeCrossSectionPerAtom( |
---|
105 | const G4ParticleDefinition* pd, |
---|
106 | G4double kinEnergy, |
---|
107 | G4double Z, |
---|
108 | G4double A, |
---|
109 | G4double cut, |
---|
110 | G4double emax) |
---|
111 | { |
---|
112 | double xs = |
---|
113 | G4KleinNishinaCompton::ComputeCrossSectionPerAtom(pd,kinEnergy, |
---|
114 | Z,A,cut,emax); |
---|
115 | G4double polzz = theBeamPolarization.p3()*theTargetPolarization.z(); |
---|
116 | if (polzz!=0) { |
---|
117 | G4double asym=ComputeAsymmetryPerAtom(kinEnergy, Z); |
---|
118 | xs*=(1.+polzz*asym); |
---|
119 | } |
---|
120 | return xs; |
---|
121 | } |
---|
122 | |
---|
123 | |
---|
124 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
125 | |
---|
126 | void G4PolarizedComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect, |
---|
127 | const G4MaterialCutsCouple*, |
---|
128 | const G4DynamicParticle* aDynamicGamma, |
---|
129 | G4double, |
---|
130 | G4double) |
---|
131 | { |
---|
132 | const G4Track * aTrack = fParticleChange->GetCurrentTrack(); |
---|
133 | G4VPhysicalVolume* aPVolume = aTrack->GetVolume(); |
---|
134 | G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume(); |
---|
135 | |
---|
136 | if (verboseLevel>=1) |
---|
137 | G4cout<<"G4PolarizedComptonModel::SampleSecondaries in " |
---|
138 | << aLVolume->GetName() <<G4endl; |
---|
139 | |
---|
140 | G4PolarizationManager * polarizationManager = G4PolarizationManager::GetInstance(); |
---|
141 | |
---|
142 | // obtain polarization of the beam |
---|
143 | theBeamPolarization = aDynamicGamma->GetPolarization(); |
---|
144 | theBeamPolarization.SetPhoton(); |
---|
145 | |
---|
146 | // obtain polarization of the media |
---|
147 | const G4bool targetIsPolarized = polarizationManager->IsPolarized(aLVolume); |
---|
148 | theTargetPolarization = polarizationManager->GetVolumePolarization(aLVolume); |
---|
149 | |
---|
150 | // if beam is linear polarized or target is transversely polarized |
---|
151 | // determine the angle to x-axis |
---|
152 | // (assumes same PRF as in the polarization definition) |
---|
153 | |
---|
154 | G4ThreeVector gamDirection0 = aDynamicGamma->GetMomentumDirection(); |
---|
155 | |
---|
156 | // transfere theTargetPolarization |
---|
157 | // into the gamma frame (problem electron is at rest) |
---|
158 | if (targetIsPolarized) |
---|
159 | theTargetPolarization.rotateUz(gamDirection0); |
---|
160 | |
---|
161 | // The scattered gamma energy is sampled according to Klein - Nishina formula. |
---|
162 | // The random number techniques of Butcher & Messel are used |
---|
163 | // (Nuc Phys 20(1960),15). |
---|
164 | // Note : Effects due to binding of atomic electrons are negliged. |
---|
165 | |
---|
166 | G4double gamEnergy0 = aDynamicGamma->GetKineticEnergy(); |
---|
167 | G4double E0_m = gamEnergy0 / electron_mass_c2 ; |
---|
168 | |
---|
169 | |
---|
170 | // |
---|
171 | // sample the energy rate of the scattered gamma |
---|
172 | // |
---|
173 | |
---|
174 | G4double epsilon, epsilonsq, onecost, sint2, greject ; |
---|
175 | |
---|
176 | G4double epsilon0 = 1./(1. + 2.*E0_m); |
---|
177 | G4double epsilon0sq = epsilon0*epsilon0; |
---|
178 | G4double alpha1 = - std::log(epsilon0); |
---|
179 | G4double alpha2 = 0.5*(1.- epsilon0sq); |
---|
180 | |
---|
181 | G4double polarization = theBeamPolarization.p3()*theTargetPolarization.p3(); |
---|
182 | do { |
---|
183 | if ( alpha1/(alpha1+alpha2) > G4UniformRand() ) { |
---|
184 | epsilon = std::exp(-alpha1*G4UniformRand()); // epsilon0**r |
---|
185 | epsilonsq = epsilon*epsilon; |
---|
186 | |
---|
187 | } else { |
---|
188 | epsilonsq = epsilon0sq + (1.- epsilon0sq)*G4UniformRand(); |
---|
189 | epsilon = std::sqrt(epsilonsq); |
---|
190 | }; |
---|
191 | |
---|
192 | onecost = (1.- epsilon)/(epsilon*E0_m); |
---|
193 | sint2 = onecost*(2.-onecost); |
---|
194 | |
---|
195 | |
---|
196 | G4double gdiced = 2.*(1./epsilon+epsilon); |
---|
197 | G4double gdist = 1./epsilon + epsilon - sint2 |
---|
198 | - polarization*(1./epsilon-epsilon)*(1.-onecost); |
---|
199 | |
---|
200 | greject = gdist/gdiced; |
---|
201 | |
---|
202 | if (greject>1) G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n" |
---|
203 | <<" costh rejection does not work properly: "<<greject<<G4endl; |
---|
204 | |
---|
205 | } while (greject < G4UniformRand()); |
---|
206 | |
---|
207 | // |
---|
208 | // scattered gamma angles. ( Z - axis along the parent gamma) |
---|
209 | // |
---|
210 | |
---|
211 | G4double cosTeta = 1. - onecost; |
---|
212 | G4double sinTeta = std::sqrt (sint2); |
---|
213 | G4double Phi; |
---|
214 | do { |
---|
215 | Phi = twopi * G4UniformRand(); |
---|
216 | G4double gdiced = 1./epsilon + epsilon - sint2 |
---|
217 | + std::abs(theBeamPolarization.p3())* |
---|
218 | ( std::abs((1./epsilon-epsilon)*cosTeta*theTargetPolarization.p3()) |
---|
219 | +(1.-epsilon)*sinTeta*(std::sqrt(sqr(theTargetPolarization.p1()) |
---|
220 | + sqr(theTargetPolarization.p2())))) |
---|
221 | +sint2*(std::sqrt(sqr(theBeamPolarization.p1()) + sqr(theBeamPolarization.p2()))); |
---|
222 | |
---|
223 | G4double gdist = 1./epsilon + epsilon - sint2 |
---|
224 | + theBeamPolarization.p3()* |
---|
225 | ((1./epsilon-epsilon)*cosTeta*theTargetPolarization.p3() |
---|
226 | +(1.-epsilon)*sinTeta*(std::cos(Phi)*theTargetPolarization.p1()+ |
---|
227 | std::sin(Phi)*theTargetPolarization.p2())) |
---|
228 | -sint2*(std::cos(2.*Phi)*theBeamPolarization.p1() |
---|
229 | +std::sin(2.*Phi)*theBeamPolarization.p2()); |
---|
230 | greject = gdist/gdiced; |
---|
231 | |
---|
232 | if (greject>1.+1.e-10 || greject<0) G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n" |
---|
233 | <<" phi rejection does not work properly: "<<greject<<G4endl; |
---|
234 | |
---|
235 | if (greject<1.e-3) { |
---|
236 | G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n" |
---|
237 | <<" phi rejection does not work properly: "<<greject<<"\n"; |
---|
238 | G4cout<<" greject="<<greject<<" phi="<<Phi<<" cost="<<cosTeta<<"\n"; |
---|
239 | G4cout<<" gdiced="<<gdiced<<" gdist="<<gdist<<"\n"; |
---|
240 | G4cout<<" eps="<<epsilon<<" 1/eps="<<1./epsilon<<"\n"; |
---|
241 | } |
---|
242 | |
---|
243 | } while (greject < G4UniformRand()); |
---|
244 | G4double dirx = sinTeta*std::cos(Phi), diry = sinTeta*std::sin(Phi), dirz = cosTeta; |
---|
245 | |
---|
246 | // |
---|
247 | // update G4VParticleChange for the scattered gamma |
---|
248 | // |
---|
249 | |
---|
250 | G4ThreeVector gamDirection1 ( dirx,diry,dirz ); |
---|
251 | gamDirection1.rotateUz(gamDirection0); |
---|
252 | G4double gamEnergy1 = epsilon*gamEnergy0; |
---|
253 | fParticleChange->SetProposedKineticEnergy(gamEnergy1); |
---|
254 | |
---|
255 | |
---|
256 | if(gamEnergy1 > lowestGammaEnergy) { |
---|
257 | fParticleChange->ProposeMomentumDirection(gamDirection1); |
---|
258 | } else { |
---|
259 | fParticleChange->ProposeTrackStatus(fStopAndKill); |
---|
260 | gamEnergy1 += fParticleChange->GetLocalEnergyDeposit(); |
---|
261 | fParticleChange->ProposeLocalEnergyDeposit(gamEnergy1); |
---|
262 | } |
---|
263 | |
---|
264 | // |
---|
265 | // kinematic of the scattered electron |
---|
266 | // |
---|
267 | |
---|
268 | G4double eKinEnergy = gamEnergy0 - gamEnergy1; |
---|
269 | G4ThreeVector eDirection = gamEnergy0*gamDirection0 - gamEnergy1*gamDirection1; |
---|
270 | eDirection = eDirection.unit(); |
---|
271 | |
---|
272 | // |
---|
273 | // calculate Stokesvector of final state photon and electron |
---|
274 | // |
---|
275 | G4ThreeVector nInteractionFrame; |
---|
276 | if((gamEnergy1 > lowestGammaEnergy) || |
---|
277 | (eKinEnergy > DBL_MIN)) { |
---|
278 | |
---|
279 | // determine interaction plane |
---|
280 | // nInteractionFrame = |
---|
281 | // G4PolarizationHelper::GetFrame(gamDirection1,eDirection); |
---|
282 | if (gamEnergy1 > lowestGammaEnergy) |
---|
283 | nInteractionFrame = G4PolarizationHelper::GetFrame(gamDirection1,gamDirection0); |
---|
284 | else |
---|
285 | nInteractionFrame = G4PolarizationHelper::GetFrame(gamDirection0, eDirection); |
---|
286 | |
---|
287 | // transfere theBeamPolarization and theTargetPolarization |
---|
288 | // into the interaction frame (note electron is in gamma frame) |
---|
289 | if (verboseLevel>=1) { |
---|
290 | G4cout << "========================================\n"; |
---|
291 | G4cout << " nInteractionFrame = " <<nInteractionFrame<<"\n"; |
---|
292 | G4cout << " GammaDirection0 = " <<gamDirection0<<"\n"; |
---|
293 | G4cout << " gammaPolarization = " <<theBeamPolarization<<"\n"; |
---|
294 | G4cout << " electronPolarization = " <<theTargetPolarization<<"\n"; |
---|
295 | } |
---|
296 | |
---|
297 | theBeamPolarization.InvRotateAz(nInteractionFrame,gamDirection0); |
---|
298 | theTargetPolarization.InvRotateAz(nInteractionFrame,gamDirection0); |
---|
299 | |
---|
300 | if (verboseLevel>=1) { |
---|
301 | G4cout << "----------------------------------------\n"; |
---|
302 | G4cout << " gammaPolarization = " <<theBeamPolarization<<"\n"; |
---|
303 | G4cout << " electronPolarization = " <<theTargetPolarization<<"\n"; |
---|
304 | G4cout << "----------------------------------------\n"; |
---|
305 | } |
---|
306 | |
---|
307 | // initialize the polarization transfer matrix |
---|
308 | crossSectionCalculator->Initialize(epsilon,E0_m,0., |
---|
309 | theBeamPolarization, |
---|
310 | theTargetPolarization,2); |
---|
311 | } |
---|
312 | |
---|
313 | // if(eKinEnergy > DBL_MIN) |
---|
314 | { |
---|
315 | // in interaction frame |
---|
316 | // calculate polarization transfer to the photon (in interaction plane) |
---|
317 | finalGammaPolarization = crossSectionCalculator->GetPol2(); |
---|
318 | if (verboseLevel>=1) G4cout << " gammaPolarization1 = " <<finalGammaPolarization<<"\n"; |
---|
319 | finalGammaPolarization.SetPhoton(); |
---|
320 | |
---|
321 | // translate polarization into particle reference frame |
---|
322 | finalGammaPolarization.RotateAz(nInteractionFrame,gamDirection1); |
---|
323 | //store polarization vector |
---|
324 | fParticleChange->ProposePolarization(finalGammaPolarization); |
---|
325 | if (finalGammaPolarization.mag() > 1.+1.e-8){ |
---|
326 | G4cout<<"ERROR in Polarizaed Compton Scattering !"<<G4endl; |
---|
327 | G4cout<<"Polarization of final photon more than 100%"<<G4endl; |
---|
328 | G4cout<<finalGammaPolarization<<" mag = "<<finalGammaPolarization.mag()<<G4endl; |
---|
329 | } |
---|
330 | if (verboseLevel>=1) { |
---|
331 | G4cout << " gammaPolarization1 = " <<finalGammaPolarization<<"\n"; |
---|
332 | G4cout << " GammaDirection1 = " <<gamDirection1<<"\n"; |
---|
333 | } |
---|
334 | } |
---|
335 | |
---|
336 | // if (ElecKineEnergy > fminimalEnergy) { |
---|
337 | { |
---|
338 | finalElectronPolarization = crossSectionCalculator->GetPol3(); |
---|
339 | if (verboseLevel>=1) |
---|
340 | G4cout << " electronPolarization1 = " <<finalElectronPolarization<<"\n"; |
---|
341 | |
---|
342 | // transfer into particle reference frame |
---|
343 | finalElectronPolarization.RotateAz(nInteractionFrame,eDirection); |
---|
344 | if (verboseLevel>=1) { |
---|
345 | G4cout << " electronPolarization1 = " <<finalElectronPolarization<<"\n"; |
---|
346 | G4cout << " ElecDirection = " <<eDirection<<"\n"; |
---|
347 | } |
---|
348 | } |
---|
349 | if (verboseLevel>=1) |
---|
350 | G4cout << "========================================\n"; |
---|
351 | |
---|
352 | |
---|
353 | if(eKinEnergy > DBL_MIN) { |
---|
354 | |
---|
355 | // create G4DynamicParticle object for the electron. |
---|
356 | G4DynamicParticle* aElectron = new G4DynamicParticle(theElectron,eDirection,eKinEnergy); |
---|
357 | //store polarization vector |
---|
358 | if (finalElectronPolarization.mag() > 1.+1.e-8){ |
---|
359 | G4cout<<"ERROR in Polarizaed Compton Scattering !"<<G4endl; |
---|
360 | G4cout<<"Polarization of final electron more than 100%"<<G4endl; |
---|
361 | G4cout<<finalElectronPolarization<<" mag = "<<finalElectronPolarization.mag()<<G4endl; |
---|
362 | } |
---|
363 | aElectron->SetPolarization(finalElectronPolarization.p1(), |
---|
364 | finalElectronPolarization.p2(), |
---|
365 | finalElectronPolarization.p3()); |
---|
366 | fvect->push_back(aElectron); |
---|
367 | } |
---|
368 | } |
---|
369 | |
---|
370 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
371 | |
---|
372 | |
---|