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