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