| 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 // File name: G4PolarizedCompton
27 //
28 // Author: Andreas Schaelicke
29 // based on code by Michel Mair
30 //
31 // Class description
32 // modified version respecting media and bea
33 // using the stokes formalism
34
35 #include "G4PolarizedCompton.hh"
36
37 #include "G4Electron.hh"
38 #include "G4EmParameters.hh"
39 #include "G4KleinNishinaCompton.hh"
40 #include "G4PhysicsTableHelper.hh"
41 #include "G4PolarizationManager.hh"
42 #include "G4PolarizedComptonModel.hh"
43 #include "G4ProductionCutsTable.hh"
44 #include "G4StokesVector.hh"
45 #include "G4SystemOfUnits.hh"
46
47 //....oooOO0OOooo........oooOO0OOooo........oo
48 G4PhysicsTable* G4PolarizedCompton::theAsymmet
49
50 G4PolarizedCompton::G4PolarizedCompton(const G
51 G4Proce
52 : G4VEmProcess(processName, type)
53 , fType(10)
54 , fBuildAsymmetryTable(true)
55 , fUseAsymmetryTable(true)
56 , fIsInitialised(false)
57 {
58 SetStartFromNullFlag(true);
59 SetBuildTableFlag(true);
60 SetSecondaryParticle(G4Electron::Electron())
61 SetProcessSubType(fComptonScattering);
62 SetMinKinEnergyPrim(1. * MeV);
63 SetSplineFlag(true);
64 fEmModel = nullptr;
65 }
66
67 //....oooOO0OOooo........oooOO0OOooo........oo
68 G4PolarizedCompton::~G4PolarizedCompton() { Cl
69
70 //....oooOO0OOooo........oooOO0OOooo........oo
71 void G4PolarizedCompton::ProcessDescription(st
72 {
73 out << "Polarized model for Compton scatteri
74
75 G4VEmProcess::ProcessDescription(out);
76 }
77
78 //....oooOO0OOooo........oooOO0OOooo........oo
79 void G4PolarizedCompton::CleanTable()
80 {
81 if(theAsymmetryTable)
82 {
83 theAsymmetryTable->clearAndDestroy();
84 delete theAsymmetryTable;
85 theAsymmetryTable = nullptr;
86 }
87 }
88
89 //....oooOO0OOooo........oooOO0OOooo........oo
90 G4bool G4PolarizedCompton::IsApplicable(const
91 {
92 return (&p == G4Gamma::Gamma());
93 }
94
95 //....oooOO0OOooo........oooOO0OOooo........oo
96 void G4PolarizedCompton::InitialiseProcess(con
97 {
98 if(!fIsInitialised)
99 {
100 fIsInitialised = true;
101 if(0 == fType)
102 {
103 if(nullptr == EmModel(0))
104 {
105 SetEmModel(new G4KleinNishinaCompton()
106 }
107 }
108 else
109 {
110 fEmModel = new G4PolarizedComptonModel()
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 }
119
120 //....oooOO0OOooo........oooOO0OOooo........oo
121 void G4PolarizedCompton::SetModel(const G4Stri
122 {
123 if(ss == "Klein-Nishina")
124 {
125 fType = 0;
126 }
127 if(ss == "Polarized-Compton")
128 {
129 fType = 10;
130 }
131 }
132
133 //....oooOO0OOooo........oooOO0OOooo........oo
134 G4double G4PolarizedCompton::GetMeanFreePath(c
135 G
136 G
137 {
138 // *** get unploarised mean free path from l
139 G4double mfp =
140 G4VEmProcess::GetMeanFreePath(aTrack, prev
141
142 if(theAsymmetryTable && fUseAsymmetryTable &
143 {
144 mfp *= ComputeSaturationFactor(aTrack);
145 }
146 if(verboseLevel >= 2)
147 {
148 G4cout << "G4PolarizedCompton::MeanFreePat
149 << G4endl;
150 }
151 return mfp;
152 }
153
154 //....oooOO0OOooo........oooOO0OOooo........oo
155 G4double G4PolarizedCompton::PostStepGetPhysic
156 const G4Track& aTrack, G4double previousStep
157 {
158 // save previous values
159 G4double nLength = theNumberOfInteractionLen
160 G4double iLength = currentInteractionLength;
161
162 // *** compute unpolarized step limit ***
163 // this changes theNumberOfInteractionLength
164 G4double x = G4VEmProcess::PostStepGetPhysic
165 aTrack, previousStepSize, condition);
166 G4double x0 = x;
167 G4double satFact = 1.0;
168
169 // *** add corrections on polarisation ***
170 if(theAsymmetryTable && fUseAsymmetryTable &
171 {
172 satFact = ComputeSaturationFact
173 G4double curLength = currentInteractionLen
174 G4double prvLength = iLength * satFact;
175 if(nLength > 0.0)
176 {
177 theNumberOfInteractionLengthLeft =
178 std::max(nLength - previousStepSize /
179 }
180 x = theNumberOfInteractionLengthLeft * cur
181 }
182 if(verboseLevel >= 2)
183 {
184 G4cout << "G4PolarizedCompton::PostStepGPI
185 << x / mm << " mm;" << G4endl
186 << " unpolarized valu
187 << x0 / mm << " mm." << G4endl;
188 }
189 return x;
190 }
191
192 //....oooOO0OOooo........oooOO0OOooo........oo
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
229 std::size_t midx = CurrentMa
230 const G4PhysicsVector* aVector = nullptr;
231 if(midx < theAsymmetryTable->size())
232 {
233 aVector = (*theAsymmetryTable)(midx);
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 }
265
266 //....oooOO0OOooo........oooOO0OOooo........oo
267 void G4PolarizedCompton::BuildPhysicsTable(con
268 {
269 // *** build (unpolarized) cross section tab
270 G4VEmProcess::BuildPhysicsTable(part);
271 if(fBuildAsymmetryTable && fEmModel)
272 {
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 }
286
287 //....oooOO0OOooo........oooOO0OOooo........oo
288 void G4PolarizedCompton::BuildAsymmetryTable(c
289 {
290 // cleanup old, initialise new table
291 CleanTable();
292 theAsymmetryTable =
293 G4PhysicsTableHelper::PreparePhysicsTable(
294
295 // Access to materials
296 const G4ProductionCutsTable* theCoupleTable
297 G4ProductionCutsTable::GetProductionCutsTa
298 G4int numOfCouples = (G4int)theCoupleTable->
299 if(!theAsymmetryTable)
300 {
301 return;
302 }
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
314 const G4MaterialCutsCouple* couple =
315 theCoupleTable->GetMaterialCutsCouple(
316 // use same parameters as for lambda
317 if(!aVector)
318 {
319 aVector = new G4PhysicsLogVector(emin,
320 bVector = aVector;
321 }
322 else
323 {
324 bVector = new G4PhysicsLogVector(*aVec
325 }
326
327 for(G4int j = 0; j <= nbins; ++j)
328 {
329 G4double energy = bVector->Energy(j);
330 G4double tasm = 0.;
331 G4double asym = ComputeAsymmetry(ene
332 bVector->PutValue(j, asym);
333 }
334 bVector->FillSecondDerivatives();
335 G4PhysicsTableHelper::SetPhysicsVector(t
336 }
337 }
338 }
339
340 //....oooOO0OOooo........oooOO0OOooo........oo
341 G4double G4PolarizedCompton::ComputeAsymmetry(
342 G4double energy, const G4MaterialCutsCouple*
343 const G4ParticleDefinition& aParticle, G4dou
344 {
345 G4double lAsymmetry = 0.0;
346 tAsymmetry = 0;
347
348 // calculate polarized cross section
349 G4ThreeVector thePolarization = G4ThreeVecto
350 fEmModel->SetTargetPolarization(thePolarizat
351 fEmModel->SetBeamPolarization(thePolarizatio
352 G4double sigma2 =
353 fEmModel->CrossSection(couple, &aParticle,
354
355 // calculate unpolarized cross section
356 thePolarization = G4ThreeVector();
357 fEmModel->SetTargetPolarization(thePolarizat
358 fEmModel->SetBeamPolarization(thePolarizatio
359 G4double sigma0 =
360 fEmModel->CrossSection(couple, &aParticle,
361
362 // determine asymmetries
363 if(sigma0 > 0.)
364 {
365 lAsymmetry = sigma2 / sigma0 - 1.;
366 }
367 return lAsymmetry;
368 }
369