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