| 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 // ------------------------------------------- 26 // -------------------------------------------------------------------
27 // 27 //
28 // Geant4 Class file 28 // Geant4 Class file
29 // 29 //
30 // File name: G4PolarizedIonisation 30 // File name: G4PolarizedIonisation
31 // 31 //
32 // Author: A.Schaelicke on base of Vlad 32 // Author: A.Schaelicke on base of Vladimir Ivanchenko code
33 33
34 #include "G4PolarizedIonisation.hh" 34 #include "G4PolarizedIonisation.hh"
35 35
36 #include "G4Electron.hh" 36 #include "G4Electron.hh"
37 #include "G4EmParameters.hh" 37 #include "G4EmParameters.hh"
38 #include "G4PhysicsTableHelper.hh" 38 #include "G4PhysicsTableHelper.hh"
39 #include "G4PolarizationHelper.hh" 39 #include "G4PolarizationHelper.hh"
40 #include "G4PolarizationManager.hh" 40 #include "G4PolarizationManager.hh"
41 #include "G4PolarizedIonisationModel.hh" 41 #include "G4PolarizedIonisationModel.hh"
42 #include "G4Positron.hh" 42 #include "G4Positron.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 #include "G4UnitsTable.hh" 46 #include "G4UnitsTable.hh"
47 #include "G4UniversalFluctuation.hh" 47 #include "G4UniversalFluctuation.hh"
48 48
49 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
50 G4PolarizedIonisation::G4PolarizedIonisation(c 50 G4PolarizedIonisation::G4PolarizedIonisation(const G4String& name)
51 : G4VEnergyLossProcess(name) 51 : G4VEnergyLossProcess(name)
52 , fAsymmetryTable(nullptr) 52 , fAsymmetryTable(nullptr)
53 , fTransverseAsymmetryTable(nullptr) 53 , fTransverseAsymmetryTable(nullptr)
54 , fIsElectron(true) 54 , fIsElectron(true)
55 , fIsInitialised(false) 55 , fIsInitialised(false)
56 { 56 {
57 verboseLevel = 0; 57 verboseLevel = 0;
58 SetProcessSubType(fIonisation); 58 SetProcessSubType(fIonisation);
59 SetSecondaryParticle(G4Electron::Electron()) 59 SetSecondaryParticle(G4Electron::Electron());
60 fFlucModel = nullptr; 60 fFlucModel = nullptr;
61 fEmModel = nullptr; 61 fEmModel = nullptr;
62 } 62 }
63 63
64 //....oooOO0OOooo........oooOO0OOooo........oo 64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
65 G4PolarizedIonisation::~G4PolarizedIonisation( 65 G4PolarizedIonisation::~G4PolarizedIonisation() { CleanTables(); }
66 66
67 //....oooOO0OOooo........oooOO0OOooo........oo 67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
68 void G4PolarizedIonisation::ProcessDescription 68 void G4PolarizedIonisation::ProcessDescription(std::ostream& out) const
69 { 69 {
70 out << "Polarized version of G4eIonisation.\ 70 out << "Polarized version of G4eIonisation.\n";
71 71
72 G4VEnergyLossProcess::ProcessDescription(out 72 G4VEnergyLossProcess::ProcessDescription(out);
73 } 73 }
74 74
75 //....oooOO0OOooo........oooOO0OOooo........oo 75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76 void G4PolarizedIonisation::CleanTables() 76 void G4PolarizedIonisation::CleanTables()
77 { 77 {
78 if(fAsymmetryTable) 78 if(fAsymmetryTable)
79 { 79 {
80 fAsymmetryTable->clearAndDestroy(); 80 fAsymmetryTable->clearAndDestroy();
81 delete fAsymmetryTable; 81 delete fAsymmetryTable;
82 fAsymmetryTable = nullptr; 82 fAsymmetryTable = nullptr;
83 } 83 }
84 if(fTransverseAsymmetryTable) 84 if(fTransverseAsymmetryTable)
85 { 85 {
86 fTransverseAsymmetryTable->clearAndDestroy 86 fTransverseAsymmetryTable->clearAndDestroy();
87 delete fTransverseAsymmetryTable; 87 delete fTransverseAsymmetryTable;
88 fTransverseAsymmetryTable = nullptr; 88 fTransverseAsymmetryTable = nullptr;
89 } 89 }
90 } 90 }
91 91
92 //....oooOO0OOooo........oooOO0OOooo........oo 92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
93 G4double G4PolarizedIonisation::MinPrimaryEner 93 G4double G4PolarizedIonisation::MinPrimaryEnergy(const G4ParticleDefinition*,
94 94 const G4Material*,
95 95 G4double cut)
96 { 96 {
97 G4double x = cut; 97 G4double x = cut;
98 if(fIsElectron) 98 if(fIsElectron)
99 { 99 {
100 x += cut; 100 x += cut;
101 } 101 }
102 return x; 102 return x;
103 } 103 }
104 104
105 //....oooOO0OOooo........oooOO0OOooo........oo 105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
106 G4bool G4PolarizedIonisation::IsApplicable(con 106 G4bool G4PolarizedIonisation::IsApplicable(const G4ParticleDefinition& p)
107 { 107 {
108 return (&p == G4Electron::Electron() || &p = 108 return (&p == G4Electron::Electron() || &p == G4Positron::Positron());
109 } 109 }
110 //....oooOO0OOooo........oooOO0OOooo........oo 110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111 void G4PolarizedIonisation::InitialiseEnergyLo 111 void G4PolarizedIonisation::InitialiseEnergyLossProcess(
112 const G4ParticleDefinition* part, const G4Pa 112 const G4ParticleDefinition* part, const G4ParticleDefinition*)
113 { 113 {
114 if(!fIsInitialised) 114 if(!fIsInitialised)
115 { 115 {
116 if(part == G4Positron::Positron()) 116 if(part == G4Positron::Positron())
117 { 117 {
118 fIsElectron = false; 118 fIsElectron = false;
119 } 119 }
120 120
121 if(!FluctModel()) 121 if(!FluctModel())
122 { 122 {
123 SetFluctModel(new G4UniversalFluctuation 123 SetFluctModel(new G4UniversalFluctuation());
124 } 124 }
125 fFlucModel = FluctModel(); 125 fFlucModel = FluctModel();
126 126
127 fEmModel = new G4PolarizedIonisationModel( 127 fEmModel = new G4PolarizedIonisationModel();
128 SetEmModel(fEmModel); 128 SetEmModel(fEmModel);
129 G4EmParameters* param = G4EmParameters::In 129 G4EmParameters* param = G4EmParameters::Instance();
130 fEmModel->SetLowEnergyLimit(param->MinKinE 130 fEmModel->SetLowEnergyLimit(param->MinKinEnergy());
131 fEmModel->SetHighEnergyLimit(param->MaxKin 131 fEmModel->SetHighEnergyLimit(param->MaxKinEnergy());
132 AddEmModel(1, fEmModel, fFlucModel); 132 AddEmModel(1, fEmModel, fFlucModel);
133 133
134 fIsInitialised = true; 134 fIsInitialised = true;
135 } 135 }
136 } 136 }
137 137
138 //....oooOO0OOooo........oooOO0OOooo........oo 138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
139 G4double G4PolarizedIonisation::GetMeanFreePat 139 G4double G4PolarizedIonisation::GetMeanFreePath(const G4Track& track,
140 140 G4double step,
141 141 G4ForceCondition* cond)
142 { 142 {
143 // *** get unploarised mean free path from l 143 // *** get unploarised mean free path from lambda table ***
144 G4double mfp = G4VEnergyLossProcess::GetMean 144 G4double mfp = G4VEnergyLossProcess::GetMeanFreePath(track, step, cond);
145 if(fAsymmetryTable && fTransverseAsymmetryTa 145 if(fAsymmetryTable && fTransverseAsymmetryTable && mfp < DBL_MAX)
146 { 146 {
147 mfp *= ComputeSaturationFactor(track); 147 mfp *= ComputeSaturationFactor(track);
148 } 148 }
149 if(verboseLevel >= 2) 149 if(verboseLevel >= 2)
150 { 150 {
151 G4cout << "G4PolarizedIonisation::MeanFree 151 G4cout << "G4PolarizedIonisation::MeanFreePath: " << mfp / mm << " mm "
152 << G4endl; 152 << G4endl;
153 } 153 }
154 return mfp; 154 return mfp;
155 } 155 }
156 156
157 //....oooOO0OOooo........oooOO0OOooo........oo 157 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
158 G4double G4PolarizedIonisation::PostStepGetPhy 158 G4double G4PolarizedIonisation::PostStepGetPhysicalInteractionLength(
159 const G4Track& track, G4double step, G4Force 159 const G4Track& track, G4double step, G4ForceCondition* cond)
160 { 160 {
161 // save previous values 161 // save previous values
162 G4double nLength = theNumberOfInteractionLen 162 G4double nLength = theNumberOfInteractionLengthLeft;
163 G4double iLength = currentInteractionLength; 163 G4double iLength = currentInteractionLength;
164 164
165 // *** get unpolarised mean free path from l 165 // *** get unpolarised mean free path from lambda table ***
166 // this changes theNumberOfInteractionLength 166 // this changes theNumberOfInteractionLengthLeft and currentInteractionLength
167 G4double x = G4VEnergyLossProcess::PostStepG 167 G4double x = G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength(
168 track, step, cond); 168 track, step, cond);
169 G4double x0 = x; 169 G4double x0 = x;
170 G4double satFact = 1.; 170 G4double satFact = 1.;
171 171
172 // *** add corrections on polarisation *** 172 // *** add corrections on polarisation ***
173 if(fAsymmetryTable && fTransverseAsymmetryTa 173 if(fAsymmetryTable && fTransverseAsymmetryTable && x < DBL_MAX)
174 { 174 {
175 satFact = ComputeSaturationFact 175 satFact = ComputeSaturationFactor(track);
176 G4double curLength = currentInteractionLen 176 G4double curLength = currentInteractionLength * satFact;
177 G4double prvLength = iLength * satFact; 177 G4double prvLength = iLength * satFact;
178 if(nLength > 0.0) 178 if(nLength > 0.0)
179 { 179 {
180 theNumberOfInteractionLengthLeft = 180 theNumberOfInteractionLengthLeft =
181 std::max(nLength - step / prvLength, 0 181 std::max(nLength - step / prvLength, 0.0);
182 } 182 }
183 x = theNumberOfInteractionLengthLeft * cur 183 x = theNumberOfInteractionLengthLeft * curLength;
184 } 184 }
185 if(verboseLevel >= 2) 185 if(verboseLevel >= 2)
186 { 186 {
187 G4cout << "G4PolarizedIonisation::PostStep 187 G4cout << "G4PolarizedIonisation::PostStepGPIL: " << std::setprecision(8)
188 << x / mm << " mm;" << G4endl 188 << x / mm << " mm;" << G4endl
189 << " unpolarized 189 << " unpolarized value: " << std::setprecision(8)
190 << x0 / mm << " mm." << G4endl; 190 << x0 / mm << " mm." << G4endl;
191 } 191 }
192 return x; 192 return x;
193 } 193 }
194 194
195 //....oooOO0OOooo........oooOO0OOooo........oo 195 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
196 G4double G4PolarizedIonisation::ComputeSaturat 196 G4double G4PolarizedIonisation::ComputeSaturationFactor(const G4Track& track)
197 { 197 {
198 const G4Material* aMaterial = track.GetMater 198 const G4Material* aMaterial = track.GetMaterial();
199 G4VPhysicalVolume* aPVolume = track.GetVolum 199 G4VPhysicalVolume* aPVolume = track.GetVolume();
200 G4LogicalVolume* aLVolume = aPVolume->GetL 200 G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
201 201
202 G4PolarizationManager* polarizationManager = 202 G4PolarizationManager* polarizationManager =
203 G4PolarizationManager::GetInstance(); 203 G4PolarizationManager::GetInstance();
204 204
205 const G4bool volumeIsPolarized = polarizatio 205 const G4bool volumeIsPolarized = polarizationManager->IsPolarized(aLVolume);
206 G4StokesVector volPolarization = 206 G4StokesVector volPolarization =
207 polarizationManager->GetVolumePolarization 207 polarizationManager->GetVolumePolarization(aLVolume);
208 208
209 G4double factor = 1.0; 209 G4double factor = 1.0;
210 210
211 if(volumeIsPolarized && !volPolarization.IsZ 211 if(volumeIsPolarized && !volPolarization.IsZero())
212 { 212 {
213 // *** get asymmetry, if target is polariz 213 // *** get asymmetry, if target is polarized ***
214 const G4DynamicParticle* aDynamicPart = tr 214 const G4DynamicParticle* aDynamicPart = track.GetDynamicParticle();
215 const G4double energy = aD 215 const G4double energy = aDynamicPart->GetKineticEnergy();
216 const G4StokesVector polarization = G4Stok 216 const G4StokesVector polarization = G4StokesVector(track.GetPolarization());
217 const G4ParticleMomentum direction0 = aDyn 217 const G4ParticleMomentum direction0 = aDynamicPart->GetMomentumDirection();
218 218
219 if(verboseLevel >= 2) 219 if(verboseLevel >= 2)
220 { 220 {
221 G4cout << "G4PolarizedIonisation::Comput 221 G4cout << "G4PolarizedIonisation::ComputeSaturationFactor: " << G4endl;
222 G4cout << " Energy(MeV) " << energy / M 222 G4cout << " Energy(MeV) " << energy / MeV << G4endl;
223 G4cout << " Direction " << direction0 223 G4cout << " Direction " << direction0 << G4endl;
224 G4cout << " Polarization " << polarizati 224 G4cout << " Polarization " << polarization << G4endl;
225 G4cout << " MaterialPol. " << volPolariz 225 G4cout << " MaterialPol. " << volPolarization << G4endl;
226 G4cout << " Phys. Volume " << aPVolume-> 226 G4cout << " Phys. Volume " << aPVolume->GetName() << G4endl;
227 G4cout << " Log. Volume " << aLVolume-> 227 G4cout << " Log. Volume " << aLVolume->GetName() << G4endl;
228 G4cout << " Material " << aMaterial 228 G4cout << " Material " << aMaterial << G4endl;
229 } 229 }
230 230
231 std::size_t midx = CurrentMa 231 std::size_t midx = CurrentMaterialCutsCoupleIndex();
232 const G4PhysicsVector* aVector = nullptr; 232 const G4PhysicsVector* aVector = nullptr;
233 const G4PhysicsVector* bVector = nullptr; 233 const G4PhysicsVector* bVector = nullptr;
234 if(midx < fAsymmetryTable->size()) 234 if(midx < fAsymmetryTable->size())
235 { 235 {
236 aVector = (*fAsymmetryTable)(midx); 236 aVector = (*fAsymmetryTable)(midx);
237 } 237 }
238 if(midx < fTransverseAsymmetryTable->size( 238 if(midx < fTransverseAsymmetryTable->size())
239 { 239 {
240 bVector = (*fTransverseAsymmetryTable)(m 240 bVector = (*fTransverseAsymmetryTable)(midx);
241 } 241 }
242 if(aVector && bVector) 242 if(aVector && bVector)
243 { 243 {
244 G4double lAsymmetry = aVector->Value(ene 244 G4double lAsymmetry = aVector->Value(energy);
245 G4double tAsymmetry = bVector->Value(ene 245 G4double tAsymmetry = bVector->Value(energy);
246 G4double polZZ = polarization.z() * 246 G4double polZZ = polarization.z() * (volPolarization * direction0);
247 G4double polXX = 247 G4double polXX =
248 polarization.x() * 248 polarization.x() *
249 (volPolarization * G4PolarizationHelpe 249 (volPolarization * G4PolarizationHelper::GetParticleFrameX(direction0));
250 G4double polYY = 250 G4double polYY =
251 polarization.y() * 251 polarization.y() *
252 (volPolarization * G4PolarizationHelpe 252 (volPolarization * G4PolarizationHelper::GetParticleFrameY(direction0));
253 253
254 factor /= (1. + polZZ * lAsymmetry + (po 254 factor /= (1. + polZZ * lAsymmetry + (polXX + polYY) * tAsymmetry);
255 255
256 if(verboseLevel >= 2) 256 if(verboseLevel >= 2)
257 { 257 {
258 G4cout << " Asymmetry: " << lAsymm 258 G4cout << " Asymmetry: " << lAsymmetry << ", " << tAsymmetry
259 << G4endl; 259 << G4endl;
260 G4cout << " PolProduct: " << polXX 260 G4cout << " PolProduct: " << polXX << ", " << polYY << ", " << polZZ
261 << G4endl; 261 << G4endl;
262 G4cout << " Factor: " << factor 262 G4cout << " Factor: " << factor << G4endl;
263 } 263 }
264 } 264 }
265 else 265 else
266 { 266 {
267 G4ExceptionDescription ed; 267 G4ExceptionDescription ed;
268 ed << "Problem with asymmetry tables: ma 268 ed << "Problem with asymmetry tables: material index " << midx
269 << " is out of range or tables are no 269 << " is out of range or tables are not filled";
270 G4Exception("G4PolarizedIonisation::Comp 270 G4Exception("G4PolarizedIonisation::ComputeSaturationFactor", "em0048",
271 JustWarning, ed, ""); 271 JustWarning, ed, "");
272 } 272 }
273 } 273 }
274 return factor; 274 return factor;
275 } 275 }
276 276
277 //....oooOO0OOooo........oooOO0OOooo........oo 277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
278 void G4PolarizedIonisation::BuildPhysicsTable( 278 void G4PolarizedIonisation::BuildPhysicsTable(const G4ParticleDefinition& part)
279 { 279 {
280 // *** build DEDX and (unpolarized) cross se 280 // *** build DEDX and (unpolarized) cross section tables
281 G4VEnergyLossProcess::BuildPhysicsTable(part 281 G4VEnergyLossProcess::BuildPhysicsTable(part);
282 G4bool master = true; 282 G4bool master = true;
283 const G4PolarizedIonisation* masterProcess = 283 const G4PolarizedIonisation* masterProcess =
284 static_cast<const G4PolarizedIonisation*>( 284 static_cast<const G4PolarizedIonisation*>(GetMasterProcess());
285 if(masterProcess && masterProcess != this) 285 if(masterProcess && masterProcess != this)
286 { 286 {
287 master = false; 287 master = false;
288 } 288 }
289 if(master) 289 if(master)
290 { 290 {
291 BuildAsymmetryTables(part); 291 BuildAsymmetryTables(part);
292 } 292 }
293 } 293 }
294 294
295 //....oooOO0OOooo........oooOO0OOooo........oo 295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
296 void G4PolarizedIonisation::BuildAsymmetryTabl 296 void G4PolarizedIonisation::BuildAsymmetryTables(
297 const G4ParticleDefinition& part) 297 const G4ParticleDefinition& part)
298 { 298 {
299 // cleanup old, initialise new table 299 // cleanup old, initialise new table
300 CleanTables(); 300 CleanTables();
301 fAsymmetryTable = G4PhysicsTableHelper::Prep 301 fAsymmetryTable = G4PhysicsTableHelper::PreparePhysicsTable(fAsymmetryTable);
302 fTransverseAsymmetryTable = 302 fTransverseAsymmetryTable =
303 G4PhysicsTableHelper::PreparePhysicsTable( 303 G4PhysicsTableHelper::PreparePhysicsTable(fTransverseAsymmetryTable);
304 304
305 const G4ProductionCutsTable* theCoupleTable 305 const G4ProductionCutsTable* theCoupleTable =
306 G4ProductionCutsTable::GetProductionCutsTa 306 G4ProductionCutsTable::GetProductionCutsTable();
307 G4int numOfCouples = (G4int)theCoupleTable-> 307 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
308 308
309 for(G4int j = 0; j < numOfCouples; ++j) 309 for(G4int j = 0; j < numOfCouples; ++j)
310 { 310 {
311 // get cut value 311 // get cut value
312 const G4MaterialCutsCouple* couple = 312 const G4MaterialCutsCouple* couple =
313 theCoupleTable->GetMaterialCutsCouple(j) 313 theCoupleTable->GetMaterialCutsCouple(j);
314 314
315 G4double cut = (*theCoupleTable->GetEnergy 315 G4double cut = (*theCoupleTable->GetEnergyCutsVector(1))[j];
316 316
317 // create physics vectors then fill it (sa 317 // create physics vectors then fill it (same parameters as lambda vector)
318 G4PhysicsVector* ptrVectorA = LambdaPhysic 318 G4PhysicsVector* ptrVectorA = LambdaPhysicsVector(couple, cut);
319 G4PhysicsVector* ptrVectorB = LambdaPhysic 319 G4PhysicsVector* ptrVectorB = LambdaPhysicsVector(couple, cut);
320 std::size_t bins = ptrVectorA-> 320 std::size_t bins = ptrVectorA->GetVectorLength();
321 321
322 for(std::size_t i = 0; i < bins; ++i) 322 for(std::size_t i = 0; i < bins; ++i)
323 { 323 {
324 G4double lowEdgeEnergy = ptrVectorA->Ene 324 G4double lowEdgeEnergy = ptrVectorA->Energy(i);
325 G4double tasm = 0.; 325 G4double tasm = 0.;
326 G4double asym = ComputeAsymmetry(lowEdge 326 G4double asym = ComputeAsymmetry(lowEdgeEnergy, couple, part, cut, tasm);
327 ptrVectorA->PutValue(i, asym); 327 ptrVectorA->PutValue(i, asym);
328 ptrVectorB->PutValue(i, tasm); 328 ptrVectorB->PutValue(i, tasm);
329 } 329 }
330 fAsymmetryTable->insertAt(j, ptrVectorA); 330 fAsymmetryTable->insertAt(j, ptrVectorA);
331 fTransverseAsymmetryTable->insertAt(j, ptr 331 fTransverseAsymmetryTable->insertAt(j, ptrVectorB);
332 } 332 }
333 } 333 }
334 334
335 //....oooOO0OOooo........oooOO0OOooo........oo 335 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
336 G4double G4PolarizedIonisation::ComputeAsymmet 336 G4double G4PolarizedIonisation::ComputeAsymmetry(
337 G4double energy, const G4MaterialCutsCouple* 337 G4double energy, const G4MaterialCutsCouple* couple,
338 const G4ParticleDefinition& aParticle, G4dou 338 const G4ParticleDefinition& aParticle, G4double cut, G4double& tAsymmetry)
339 { 339 {
340 G4double lAsymmetry = 0.0; 340 G4double lAsymmetry = 0.0;
341 tAsymmetry = 0.0; 341 tAsymmetry = 0.0;
342 if(fIsElectron) 342 if(fIsElectron)
343 { 343 {
344 lAsymmetry = tAsymmetry = -1.0; 344 lAsymmetry = tAsymmetry = -1.0;
345 } 345 }
346 346
347 // calculate polarized cross section 347 // calculate polarized cross section
348 G4ThreeVector targetPolarization = G4ThreeVe 348 G4ThreeVector targetPolarization = G4ThreeVector(0., 0., 1.);
349 fEmModel->SetTargetPolarization(targetPolari 349 fEmModel->SetTargetPolarization(targetPolarization);
350 fEmModel->SetBeamPolarization(targetPolariza 350 fEmModel->SetBeamPolarization(targetPolarization);
351 G4double sigma2 = 351 G4double sigma2 =
352 fEmModel->CrossSection(couple, &aParticle, 352 fEmModel->CrossSection(couple, &aParticle, energy, cut, energy);
353 353
354 // calculate transversely polarized cross se 354 // calculate transversely polarized cross section
355 targetPolarization = G4ThreeVector(1., 0., 0 355 targetPolarization = G4ThreeVector(1., 0., 0.);
356 fEmModel->SetTargetPolarization(targetPolari 356 fEmModel->SetTargetPolarization(targetPolarization);
357 fEmModel->SetBeamPolarization(targetPolariza 357 fEmModel->SetBeamPolarization(targetPolarization);
358 G4double sigma3 = 358 G4double sigma3 =
359 fEmModel->CrossSection(couple, &aParticle, 359 fEmModel->CrossSection(couple, &aParticle, energy, cut, energy);
360 360
361 // calculate unpolarized cross section 361 // calculate unpolarized cross section
362 targetPolarization = G4ThreeVector(); 362 targetPolarization = G4ThreeVector();
363 fEmModel->SetTargetPolarization(targetPolari 363 fEmModel->SetTargetPolarization(targetPolarization);
364 fEmModel->SetBeamPolarization(targetPolariza 364 fEmModel->SetBeamPolarization(targetPolarization);
365 G4double sigma0 = 365 G4double sigma0 =
366 fEmModel->CrossSection(couple, &aParticle, 366 fEmModel->CrossSection(couple, &aParticle, energy, cut, energy);
367 // determine asymmetries 367 // determine asymmetries
368 if(sigma0 > 0.) 368 if(sigma0 > 0.)
369 { 369 {
370 lAsymmetry = sigma2 / sigma0 - 1.; 370 lAsymmetry = sigma2 / sigma0 - 1.;
371 tAsymmetry = sigma3 / sigma0 - 1.; 371 tAsymmetry = sigma3 / sigma0 - 1.;
372 } 372 }
373 if(std::fabs(lAsymmetry) > 1.) 373 if(std::fabs(lAsymmetry) > 1.)
374 { 374 {
375 G4ExceptionDescription ed; 375 G4ExceptionDescription ed;
376 ed << "G4PolarizedIonisation::ComputeAsymm 376 ed << "G4PolarizedIonisation::ComputeAsymmetry : E(MeV)= " << energy
377 << " lAsymmetry= " << lAsymmetry << " ( 377 << " lAsymmetry= " << lAsymmetry << " (" << std::fabs(lAsymmetry) - 1.
378 << ")"; 378 << ")";
379 G4Exception("G4PolarizedIonisation::Comput 379 G4Exception("G4PolarizedIonisation::ComputeAsymmetry", "pol002",
380 JustWarning, ed); 380 JustWarning, ed);
381 } 381 }
382 if(std::fabs(tAsymmetry) > 1.) 382 if(std::fabs(tAsymmetry) > 1.)
383 { 383 {
384 G4ExceptionDescription ed; 384 G4ExceptionDescription ed;
385 ed << "G4PolarizedIonisation::ComputeAsymm 385 ed << "G4PolarizedIonisation::ComputeAsymmetry : E(MeV)= " << energy
386 << " tAsymmetry= " << tAsymmetry << " ( 386 << " tAsymmetry= " << tAsymmetry << " (" << std::fabs(tAsymmetry) - 1.
387 << ")"; 387 << ")";
388 G4Exception("G4PolarizedIonisation::Comput 388 G4Exception("G4PolarizedIonisation::ComputeAsymmetry", "pol003",
389 JustWarning, ed); 389 JustWarning, ed);
390 } 390 }
391 return lAsymmetry; 391 return lAsymmetry;
392 } 392 }
393 393