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