| 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 // File name: G4PolarizedCompton 26 // File name: G4PolarizedCompton
27 // 27 //
28 // Author: Andreas Schaelicke 28 // Author: Andreas Schaelicke
29 // based on code by Michel Mair 29 // based on code by Michel Maire / Vladimir IVANTCHENKO
30 // 30 //
31 // Class description 31 // Class description
32 // modified version respecting media and bea 32 // modified version respecting media and beam polarization
33 // using the stokes formalism 33 // using the stokes formalism
34 34
35 #include "G4PolarizedCompton.hh" 35 #include "G4PolarizedCompton.hh"
36 36
37 #include "G4Electron.hh" 37 #include "G4Electron.hh"
38 #include "G4EmParameters.hh" 38 #include "G4EmParameters.hh"
39 #include "G4KleinNishinaCompton.hh" 39 #include "G4KleinNishinaCompton.hh"
40 #include "G4PhysicsTableHelper.hh" 40 #include "G4PhysicsTableHelper.hh"
41 #include "G4PolarizationManager.hh" 41 #include "G4PolarizationManager.hh"
42 #include "G4PolarizedComptonModel.hh" 42 #include "G4PolarizedComptonModel.hh"
43 #include "G4ProductionCutsTable.hh" 43 #include "G4ProductionCutsTable.hh"
44 #include "G4StokesVector.hh" 44 #include "G4StokesVector.hh"
45 #include "G4SystemOfUnits.hh" 45 #include "G4SystemOfUnits.hh"
46 46
47 //....oooOO0OOooo........oooOO0OOooo........oo 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
48 G4PhysicsTable* G4PolarizedCompton::theAsymmet 48 G4PhysicsTable* G4PolarizedCompton::theAsymmetryTable = nullptr;
49 49
50 G4PolarizedCompton::G4PolarizedCompton(const G 50 G4PolarizedCompton::G4PolarizedCompton(const G4String& processName,
51 G4Proce 51 G4ProcessType type)
52 : G4VEmProcess(processName, type) 52 : G4VEmProcess(processName, type)
53 , fType(10) 53 , fType(10)
54 , fBuildAsymmetryTable(true) 54 , fBuildAsymmetryTable(true)
55 , fUseAsymmetryTable(true) 55 , fUseAsymmetryTable(true)
56 , fIsInitialised(false) 56 , fIsInitialised(false)
57 { 57 {
58 SetStartFromNullFlag(true); 58 SetStartFromNullFlag(true);
59 SetBuildTableFlag(true); 59 SetBuildTableFlag(true);
60 SetSecondaryParticle(G4Electron::Electron()) 60 SetSecondaryParticle(G4Electron::Electron());
61 SetProcessSubType(fComptonScattering); 61 SetProcessSubType(fComptonScattering);
62 SetMinKinEnergyPrim(1. * MeV); 62 SetMinKinEnergyPrim(1. * MeV);
63 SetSplineFlag(true); 63 SetSplineFlag(true);
64 fEmModel = nullptr; 64 fEmModel = nullptr;
65 } 65 }
66 66
67 //....oooOO0OOooo........oooOO0OOooo........oo 67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
68 G4PolarizedCompton::~G4PolarizedCompton() { Cl 68 G4PolarizedCompton::~G4PolarizedCompton() { CleanTable(); }
69 69
70 //....oooOO0OOooo........oooOO0OOooo........oo 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71 void G4PolarizedCompton::ProcessDescription(st 71 void G4PolarizedCompton::ProcessDescription(std::ostream& out) const
72 { 72 {
73 out << "Polarized model for Compton scatteri 73 out << "Polarized model for Compton scattering.\n";
74 74
75 G4VEmProcess::ProcessDescription(out); 75 G4VEmProcess::ProcessDescription(out);
76 } 76 }
77 77
78 //....oooOO0OOooo........oooOO0OOooo........oo 78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
79 void G4PolarizedCompton::CleanTable() 79 void G4PolarizedCompton::CleanTable()
80 { 80 {
81 if(theAsymmetryTable) 81 if(theAsymmetryTable)
82 { 82 {
83 theAsymmetryTable->clearAndDestroy(); 83 theAsymmetryTable->clearAndDestroy();
84 delete theAsymmetryTable; 84 delete theAsymmetryTable;
85 theAsymmetryTable = nullptr; 85 theAsymmetryTable = nullptr;
86 } 86 }
87 } 87 }
88 88
89 //....oooOO0OOooo........oooOO0OOooo........oo 89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
90 G4bool G4PolarizedCompton::IsApplicable(const 90 G4bool G4PolarizedCompton::IsApplicable(const G4ParticleDefinition& p)
91 { 91 {
92 return (&p == G4Gamma::Gamma()); 92 return (&p == G4Gamma::Gamma());
93 } 93 }
94 94
95 //....oooOO0OOooo........oooOO0OOooo........oo 95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
96 void G4PolarizedCompton::InitialiseProcess(con 96 void G4PolarizedCompton::InitialiseProcess(const G4ParticleDefinition*)
97 { 97 {
98 if(!fIsInitialised) 98 if(!fIsInitialised)
99 { 99 {
100 fIsInitialised = true; 100 fIsInitialised = true;
101 if(0 == fType) 101 if(0 == fType)
102 { 102 {
103 if(nullptr == EmModel(0)) 103 if(nullptr == EmModel(0))
104 { 104 {
105 SetEmModel(new G4KleinNishinaCompton() 105 SetEmModel(new G4KleinNishinaCompton());
106 } 106 }
107 } 107 }
108 else 108 else
109 { 109 {
110 fEmModel = new G4PolarizedComptonModel() 110 fEmModel = new G4PolarizedComptonModel();
111 SetEmModel(fEmModel); 111 SetEmModel(fEmModel);
112 } 112 }
113 G4EmParameters* param = G4EmParameters::In 113 G4EmParameters* param = G4EmParameters::Instance();
114 EmModel(0)->SetLowEnergyLimit(param->MinKi 114 EmModel(0)->SetLowEnergyLimit(param->MinKinEnergy());
115 EmModel(0)->SetHighEnergyLimit(param->MaxK 115 EmModel(0)->SetHighEnergyLimit(param->MaxKinEnergy());
116 AddEmModel(1, EmModel(0)); 116 AddEmModel(1, EmModel(0));
117 } 117 }
118 } 118 }
119 119
120 //....oooOO0OOooo........oooOO0OOooo........oo 120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
121 void G4PolarizedCompton::SetModel(const G4Stri 121 void G4PolarizedCompton::SetModel(const G4String& ss)
122 { 122 {
123 if(ss == "Klein-Nishina") 123 if(ss == "Klein-Nishina")
124 { 124 {
125 fType = 0; 125 fType = 0;
126 } 126 }
127 if(ss == "Polarized-Compton") 127 if(ss == "Polarized-Compton")
128 { 128 {
129 fType = 10; 129 fType = 10;
130 } 130 }
131 } 131 }
132 132
133 //....oooOO0OOooo........oooOO0OOooo........oo 133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
134 G4double G4PolarizedCompton::GetMeanFreePath(c 134 G4double G4PolarizedCompton::GetMeanFreePath(const G4Track& aTrack,
135 G 135 G4double previousStepSize,
136 G 136 G4ForceCondition* condition)
137 { 137 {
138 // *** get unploarised mean free path from l 138 // *** get unploarised mean free path from lambda table ***
139 G4double mfp = 139 G4double mfp =
140 G4VEmProcess::GetMeanFreePath(aTrack, prev 140 G4VEmProcess::GetMeanFreePath(aTrack, previousStepSize, condition);
141 141
142 if(theAsymmetryTable && fUseAsymmetryTable & 142 if(theAsymmetryTable && fUseAsymmetryTable && mfp < DBL_MAX)
143 { 143 {
144 mfp *= ComputeSaturationFactor(aTrack); 144 mfp *= ComputeSaturationFactor(aTrack);
145 } 145 }
146 if(verboseLevel >= 2) 146 if(verboseLevel >= 2)
147 { 147 {
148 G4cout << "G4PolarizedCompton::MeanFreePat 148 G4cout << "G4PolarizedCompton::MeanFreePath: " << mfp / mm << " mm "
149 << G4endl; 149 << G4endl;
150 } 150 }
151 return mfp; 151 return mfp;
152 } 152 }
153 153
154 //....oooOO0OOooo........oooOO0OOooo........oo 154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
155 G4double G4PolarizedCompton::PostStepGetPhysic 155 G4double G4PolarizedCompton::PostStepGetPhysicalInteractionLength(
156 const G4Track& aTrack, G4double previousStep 156 const G4Track& aTrack, G4double previousStepSize, G4ForceCondition* condition)
157 { 157 {
158 // save previous values 158 // save previous values
159 G4double nLength = theNumberOfInteractionLen 159 G4double nLength = theNumberOfInteractionLengthLeft;
160 G4double iLength = currentInteractionLength; 160 G4double iLength = currentInteractionLength;
161 161
162 // *** compute unpolarized step limit *** 162 // *** compute unpolarized step limit ***
163 // this changes theNumberOfInteractionLength 163 // this changes theNumberOfInteractionLengthLeft and currentInteractionLength
164 G4double x = G4VEmProcess::PostStepGetPhysic 164 G4double x = G4VEmProcess::PostStepGetPhysicalInteractionLength(
165 aTrack, previousStepSize, condition); 165 aTrack, previousStepSize, condition);
166 G4double x0 = x; 166 G4double x0 = x;
167 G4double satFact = 1.0; 167 G4double satFact = 1.0;
168 168
169 // *** add corrections on polarisation *** 169 // *** add corrections on polarisation ***
170 if(theAsymmetryTable && fUseAsymmetryTable & 170 if(theAsymmetryTable && fUseAsymmetryTable && x < DBL_MAX)
171 { 171 {
172 satFact = ComputeSaturationFact 172 satFact = ComputeSaturationFactor(aTrack);
173 G4double curLength = currentInteractionLen 173 G4double curLength = currentInteractionLength * satFact;
174 G4double prvLength = iLength * satFact; 174 G4double prvLength = iLength * satFact;
175 if(nLength > 0.0) 175 if(nLength > 0.0)
176 { 176 {
177 theNumberOfInteractionLengthLeft = 177 theNumberOfInteractionLengthLeft =
178 std::max(nLength - previousStepSize / 178 std::max(nLength - previousStepSize / prvLength, 0.0);
179 } 179 }
180 x = theNumberOfInteractionLengthLeft * cur 180 x = theNumberOfInteractionLengthLeft * curLength;
181 } 181 }
182 if(verboseLevel >= 2) 182 if(verboseLevel >= 2)
183 { 183 {
184 G4cout << "G4PolarizedCompton::PostStepGPI 184 G4cout << "G4PolarizedCompton::PostStepGPIL: " << std::setprecision(8)
185 << x / mm << " mm;" << G4endl 185 << x / mm << " mm;" << G4endl
186 << " unpolarized valu 186 << " unpolarized value: " << std::setprecision(8)
187 << x0 / mm << " mm." << G4endl; 187 << x0 / mm << " mm." << G4endl;
188 } 188 }
189 return x; 189 return x;
190 } 190 }
191 191
192 //....oooOO0OOooo........oooOO0OOooo........oo 192 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
193 G4double G4PolarizedCompton::ComputeSaturation 193 G4double G4PolarizedCompton::ComputeSaturationFactor(const G4Track& aTrack)
194 { 194 {
195 G4double factor = 1.0; 195 G4double factor = 1.0;
196 196
197 // *** get asymmetry, if target is polarized 197 // *** get asymmetry, if target is polarized ***
198 const G4DynamicParticle* aDynamicGamma = aTr 198 const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
199 const G4double GammaEnergy = aDy 199 const G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
200 const G4StokesVector GammaPolarization = 200 const G4StokesVector GammaPolarization =
201 G4StokesVector(aTrack.GetPolarization()); 201 G4StokesVector(aTrack.GetPolarization());
202 const G4ParticleMomentum GammaDirection0 = 202 const G4ParticleMomentum GammaDirection0 =
203 aDynamicGamma->GetMomentumDirection(); 203 aDynamicGamma->GetMomentumDirection();
204 204
205 const G4Material* aMaterial = aTrack.GetMate 205 const G4Material* aMaterial = aTrack.GetMaterial();
206 G4VPhysicalVolume* aPVolume = aTrack.GetVolu 206 G4VPhysicalVolume* aPVolume = aTrack.GetVolume();
207 G4LogicalVolume* aLVolume = aPVolume->GetL 207 G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
208 208
209 G4PolarizationManager* polarizationManager = 209 G4PolarizationManager* polarizationManager =
210 G4PolarizationManager::GetInstance(); 210 G4PolarizationManager::GetInstance();
211 211
212 const G4bool VolumeIsPolarized = polarizatio 212 const G4bool VolumeIsPolarized = polarizationManager->IsPolarized(aLVolume);
213 G4StokesVector ElectronPolarization = 213 G4StokesVector ElectronPolarization =
214 polarizationManager->GetVolumePolarization 214 polarizationManager->GetVolumePolarization(aLVolume);
215 215
216 if(VolumeIsPolarized) 216 if(VolumeIsPolarized)
217 { 217 {
218 if(verboseLevel >= 2) 218 if(verboseLevel >= 2)
219 { 219 {
220 G4cout << "G4PolarizedCompton::ComputeSa 220 G4cout << "G4PolarizedCompton::ComputeSaturationFactor: " << G4endl;
221 G4cout << " Mom " << GammaDirection0 << 221 G4cout << " Mom " << GammaDirection0 << G4endl;
222 G4cout << " Polarization " << GammaPolar 222 G4cout << " Polarization " << GammaPolarization << G4endl;
223 G4cout << " MaterialPol. " << ElectronPo 223 G4cout << " MaterialPol. " << ElectronPolarization << G4endl;
224 G4cout << " Phys. Volume " << aPVolume-> 224 G4cout << " Phys. Volume " << aPVolume->GetName() << G4endl;
225 G4cout << " Log. Volume " << aLVolume-> 225 G4cout << " Log. Volume " << aLVolume->GetName() << G4endl;
226 G4cout << " Material " << aMaterial 226 G4cout << " Material " << aMaterial << G4endl;
227 } 227 }
228 228
229 std::size_t midx = CurrentMa 229 std::size_t midx = CurrentMaterialCutsCoupleIndex();
230 const G4PhysicsVector* aVector = nullptr; 230 const G4PhysicsVector* aVector = nullptr;
231 if(midx < theAsymmetryTable->size()) 231 if(midx < theAsymmetryTable->size())
232 { 232 {
233 aVector = (*theAsymmetryTable)(midx); 233 aVector = (*theAsymmetryTable)(midx);
234 } 234 }
235 if(aVector) 235 if(aVector)
236 { 236 {
237 G4double asymmetry = aVector->Value(Gamm 237 G4double asymmetry = aVector->Value(GammaEnergy);
238 238
239 // we have to determine angle between p 239 // we have to determine angle between particle motion
240 // and target polarisation here 240 // and target polarisation here
241 // circ pol * Vec(ElectronPol)*Vec( 241 // circ pol * Vec(ElectronPol)*Vec(PhotonMomentum)
242 // both vectors in global reference fra 242 // both vectors in global reference frame
243 243
244 G4double pol = ElectronPolarizati 244 G4double pol = ElectronPolarization * GammaDirection0;
245 G4double polProduct = GammaPolarization. 245 G4double polProduct = GammaPolarization.p3() * pol;
246 factor /= (1. + polProduct * asymmetry); 246 factor /= (1. + polProduct * asymmetry);
247 if(verboseLevel >= 2) 247 if(verboseLevel >= 2)
248 { 248 {
249 G4cout << " Asymmetry: " << asymme 249 G4cout << " Asymmetry: " << asymmetry << G4endl;
250 G4cout << " PolProduct: " << polPro 250 G4cout << " PolProduct: " << polProduct << G4endl;
251 G4cout << " Factor: " << factor 251 G4cout << " Factor: " << factor << G4endl;
252 } 252 }
253 } 253 }
254 else 254 else
255 { 255 {
256 G4ExceptionDescription ed; 256 G4ExceptionDescription ed;
257 ed << "Problem with asymmetry table: mat 257 ed << "Problem with asymmetry table: material index " << midx
258 << " is out of range or the table is 258 << " is out of range or the table is not filled";
259 G4Exception("G4PolarizedComptonModel::Co 259 G4Exception("G4PolarizedComptonModel::ComputeSaturationFactor", "em0048",
260 JustWarning, ed, ""); 260 JustWarning, ed, "");
261 } 261 }
262 } 262 }
263 return factor; 263 return factor;
264 } 264 }
265 265
266 //....oooOO0OOooo........oooOO0OOooo........oo 266 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
267 void G4PolarizedCompton::BuildPhysicsTable(con 267 void G4PolarizedCompton::BuildPhysicsTable(const G4ParticleDefinition& part)
268 { 268 {
269 // *** build (unpolarized) cross section tab 269 // *** build (unpolarized) cross section tables (Lambda)
270 G4VEmProcess::BuildPhysicsTable(part); 270 G4VEmProcess::BuildPhysicsTable(part);
271 if(fBuildAsymmetryTable && fEmModel) 271 if(fBuildAsymmetryTable && fEmModel)
272 { 272 {
273 G4bool isMaster = true; 273 G4bool isMaster = true;
274 const G4PolarizedCompton* masterProcess = 274 const G4PolarizedCompton* masterProcess =
275 static_cast<const G4PolarizedCompton*>(G 275 static_cast<const G4PolarizedCompton*>(GetMasterProcess());
276 if(masterProcess && masterProcess != this) 276 if(masterProcess && masterProcess != this)
277 { 277 {
278 isMaster = false; 278 isMaster = false;
279 } 279 }
280 if(isMaster) 280 if(isMaster)
281 { 281 {
282 BuildAsymmetryTable(part); 282 BuildAsymmetryTable(part);
283 } 283 }
284 } 284 }
285 } 285 }
286 286
287 //....oooOO0OOooo........oooOO0OOooo........oo 287 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
288 void G4PolarizedCompton::BuildAsymmetryTable(c 288 void G4PolarizedCompton::BuildAsymmetryTable(const G4ParticleDefinition& part)
289 { 289 {
290 // cleanup old, initialise new table 290 // cleanup old, initialise new table
291 CleanTable(); 291 CleanTable();
292 theAsymmetryTable = 292 theAsymmetryTable =
293 G4PhysicsTableHelper::PreparePhysicsTable( 293 G4PhysicsTableHelper::PreparePhysicsTable(theAsymmetryTable);
294 294
295 // Access to materials 295 // Access to materials
296 const G4ProductionCutsTable* theCoupleTable 296 const G4ProductionCutsTable* theCoupleTable =
297 G4ProductionCutsTable::GetProductionCutsTa 297 G4ProductionCutsTable::GetProductionCutsTable();
298 G4int numOfCouples = (G4int)theCoupleTable-> 298 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
299 if(!theAsymmetryTable) 299 if(!theAsymmetryTable)
300 { 300 {
301 return; 301 return;
302 } 302 }
303 G4int nbins = LambdaBinning( 303 G4int nbins = LambdaBinning();
304 G4double emin = MinKinEnergy() 304 G4double emin = MinKinEnergy();
305 G4double emax = MaxKinEnergy() 305 G4double emax = MaxKinEnergy();
306 G4PhysicsLogVector* aVector = nullptr; 306 G4PhysicsLogVector* aVector = nullptr;
307 G4PhysicsLogVector* bVector = nullptr; 307 G4PhysicsLogVector* bVector = nullptr;
308 308
309 for(G4int i = 0; i < numOfCouples; ++i) 309 for(G4int i = 0; i < numOfCouples; ++i)
310 { 310 {
311 if(theAsymmetryTable->GetFlag(i)) 311 if(theAsymmetryTable->GetFlag(i))
312 { 312 {
313 // create physics vector and fill it 313 // create physics vector and fill it
314 const G4MaterialCutsCouple* couple = 314 const G4MaterialCutsCouple* couple =
315 theCoupleTable->GetMaterialCutsCouple( 315 theCoupleTable->GetMaterialCutsCouple(i);
316 // use same parameters as for lambda 316 // use same parameters as for lambda
317 if(!aVector) 317 if(!aVector)
318 { 318 {
319 aVector = new G4PhysicsLogVector(emin, 319 aVector = new G4PhysicsLogVector(emin, emax, nbins, true);
320 bVector = aVector; 320 bVector = aVector;
321 } 321 }
322 else 322 else
323 { 323 {
324 bVector = new G4PhysicsLogVector(*aVec 324 bVector = new G4PhysicsLogVector(*aVector);
325 } 325 }
326 326
327 for(G4int j = 0; j <= nbins; ++j) 327 for(G4int j = 0; j <= nbins; ++j)
328 { 328 {
329 G4double energy = bVector->Energy(j); 329 G4double energy = bVector->Energy(j);
330 G4double tasm = 0.; 330 G4double tasm = 0.;
331 G4double asym = ComputeAsymmetry(ene 331 G4double asym = ComputeAsymmetry(energy, couple, part, 0., tasm);
332 bVector->PutValue(j, asym); 332 bVector->PutValue(j, asym);
333 } 333 }
334 bVector->FillSecondDerivatives(); 334 bVector->FillSecondDerivatives();
335 G4PhysicsTableHelper::SetPhysicsVector(t 335 G4PhysicsTableHelper::SetPhysicsVector(theAsymmetryTable, i, bVector);
336 } 336 }
337 } 337 }
338 } 338 }
339 339
340 //....oooOO0OOooo........oooOO0OOooo........oo 340 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
341 G4double G4PolarizedCompton::ComputeAsymmetry( 341 G4double G4PolarizedCompton::ComputeAsymmetry(
342 G4double energy, const G4MaterialCutsCouple* 342 G4double energy, const G4MaterialCutsCouple* couple,
343 const G4ParticleDefinition& aParticle, G4dou 343 const G4ParticleDefinition& aParticle, G4double cut, G4double& tAsymmetry)
344 { 344 {
345 G4double lAsymmetry = 0.0; 345 G4double lAsymmetry = 0.0;
346 tAsymmetry = 0; 346 tAsymmetry = 0;
347 347
348 // calculate polarized cross section 348 // calculate polarized cross section
349 G4ThreeVector thePolarization = G4ThreeVecto 349 G4ThreeVector thePolarization = G4ThreeVector(0., 0., 1.);
350 fEmModel->SetTargetPolarization(thePolarizat 350 fEmModel->SetTargetPolarization(thePolarization);
351 fEmModel->SetBeamPolarization(thePolarizatio 351 fEmModel->SetBeamPolarization(thePolarization);
352 G4double sigma2 = 352 G4double sigma2 =
353 fEmModel->CrossSection(couple, &aParticle, 353 fEmModel->CrossSection(couple, &aParticle, energy, cut, energy);
354 354
355 // calculate unpolarized cross section 355 // calculate unpolarized cross section
356 thePolarization = G4ThreeVector(); 356 thePolarization = G4ThreeVector();
357 fEmModel->SetTargetPolarization(thePolarizat 357 fEmModel->SetTargetPolarization(thePolarization);
358 fEmModel->SetBeamPolarization(thePolarizatio 358 fEmModel->SetBeamPolarization(thePolarization);
359 G4double sigma0 = 359 G4double sigma0 =
360 fEmModel->CrossSection(couple, &aParticle, 360 fEmModel->CrossSection(couple, &aParticle, energy, cut, energy);
361 361
362 // determine asymmetries 362 // determine asymmetries
363 if(sigma0 > 0.) 363 if(sigma0 > 0.)
364 { 364 {
365 lAsymmetry = sigma2 / sigma0 - 1.; 365 lAsymmetry = sigma2 / sigma0 - 1.;
366 } 366 }
367 return lAsymmetry; 367 return lAsymmetry;
368 } 368 }
369 369