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