ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PolarizedCompton.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4PolarizedCompton.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 //
28 //
29 // File name: G4PolarizedCompton
30 //
31 // Author: Andreas Schaelicke
32 // based on code by Michel Maire / Vladimir IVANTCHENKO
33 // Class description
34 //
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 // -----------------------------------------------------------------------------
53 
54 
55 #include "G4PolarizedCompton.hh"
56 #include "G4SystemOfUnits.hh"
57 #include "G4Electron.hh"
58 
59 #include "G4StokesVector.hh"
60 #include "G4PolarizationManager.hh"
61 #include "G4PolarizedComptonModel.hh"
62 #include "G4ProductionCutsTable.hh"
63 #include "G4PhysicsTableHelper.hh"
64 #include "G4KleinNishinaCompton.hh"
65 #include "G4PolarizedComptonModel.hh"
66 #include "G4EmParameters.hh"
67 
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
69 
70 G4PhysicsTable* G4PolarizedCompton::theAsymmetryTable = nullptr;
71 
72 G4PolarizedCompton::G4PolarizedCompton(const G4String& processName,
73  G4ProcessType type):
74  G4VEmProcess (processName, type),
75  buildAsymmetryTable(true),
76  useAsymmetryTable(true),
77  isInitialised(false),
78  mType(10),
79  targetPolarization(0.0,0.0,0.0)
80 {
81  SetStartFromNullFlag(true);
82  SetBuildTableFlag(true);
83  SetSecondaryParticle(G4Electron::Electron());
84  SetProcessSubType(fComptonScattering);
85  SetMinKinEnergyPrim(1*MeV);
86  SetSplineFlag(true);
87  emModel = nullptr;
88 }
89 
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
91 
92 G4PolarizedCompton::~G4PolarizedCompton()
93 {
94  CleanTable();
95 }
96 
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
98 
99 void G4PolarizedCompton::CleanTable()
100 {
101  if( theAsymmetryTable) {
102  theAsymmetryTable->clearAndDestroy();
103  delete theAsymmetryTable;
104  theAsymmetryTable = nullptr;
105  }
106 }
107 
108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
109 
110 G4bool G4PolarizedCompton::IsApplicable(const G4ParticleDefinition& p)
111 {
112  return (&p == G4Gamma::Gamma());
113 }
114 
115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
116 
117 void G4PolarizedCompton::InitialiseProcess(const G4ParticleDefinition*)
118 {
119  if(!isInitialised) {
120  isInitialised = true;
121  if(0 == mType) {
122  if(!EmModel(0)) { SetEmModel(new G4KleinNishinaCompton()); }
123  } else {
124  emModel = new G4PolarizedComptonModel();
125  SetEmModel(emModel);
126  }
127  G4EmParameters* param = G4EmParameters::Instance();
128  EmModel(0)->SetLowEnergyLimit(param->MinKinEnergy());
129  EmModel(0)->SetHighEnergyLimit(param->MaxKinEnergy());
130  AddEmModel(1, EmModel(0));
131  }
132 }
133 
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 
146 void G4PolarizedCompton::SetModel(const G4String& ss)
147 {
148  if(ss == "Klein-Nishina") { mType = 0; }
149  if(ss == "Polarized-Compton") { mType = 10; }
150 }
151 
152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
153 
154 G4double G4PolarizedCompton::GetMeanFreePath(const G4Track& aTrack,
155  G4double previousStepSize,
156  G4ForceCondition* condition)
157 {
158  // *** get unploarised mean free path from lambda table ***
159  G4double mfp = G4VEmProcess::GetMeanFreePath(aTrack, previousStepSize, condition);
160 
161  if (theAsymmetryTable && useAsymmetryTable && mfp < DBL_MAX) {
162  mfp *= ComputeSaturationFactor(aTrack);
163  }
164  if (verboseLevel>=2) {
165  G4cout << "G4PolarizedCompton::MeanFreePath: " << mfp / mm << " mm " << G4endl;
166  }
167  return mfp;
168 }
169 
170 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
171 
172 G4double G4PolarizedCompton::PostStepGetPhysicalInteractionLength(
173  const G4Track& aTrack,
174  G4double previousStepSize,
175  G4ForceCondition* condition)
176 {
177  // save previous values
178  G4double nLength = theNumberOfInteractionLengthLeft;
179  G4double iLength = currentInteractionLength;
180 
181  // *** compute unpolarized step limit ***
182  // this changes theNumberOfInteractionLengthLeft and currentInteractionLength
183  G4double x = G4VEmProcess::PostStepGetPhysicalInteractionLength(aTrack,
184  previousStepSize,
185  condition);
186  G4double x0 = x;
187  G4double satFact = 1.0;
188 
189  // *** add corrections on polarisation ***
190  if (theAsymmetryTable && useAsymmetryTable && x < DBL_MAX) {
191  satFact = ComputeSaturationFactor(aTrack);
192  G4double curLength = currentInteractionLength*satFact;
193  G4double prvLength = iLength*satFact;
194  if(nLength > 0.0) {
195  theNumberOfInteractionLengthLeft =
196  std::max(nLength - previousStepSize/prvLength, 0.0);
197  }
198  x = theNumberOfInteractionLengthLeft * curLength;
199  }
200  if (verboseLevel>=2) {
201  G4cout << "G4PolarizedCompton::PostStepGPIL: "
202  << std::setprecision(8) << x/mm << " mm;" << G4endl
203  << " unpolarized value: "
204  << std::setprecision(8) << x0/mm << " mm." << G4endl;
205  }
206  return x;
207 }
208 
209 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
210 
211 G4double G4PolarizedCompton::ComputeSaturationFactor(const G4Track& aTrack)
212 {
213  G4double factor = 1.0;
214 
215  // *** get asymmetry, if target is polarized ***
216  const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
217  const G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
218  const G4StokesVector GammaPolarization = aTrack.GetPolarization();
219  const G4ParticleMomentum GammaDirection0 = aDynamicGamma->GetMomentumDirection();
220 
221  G4Material* aMaterial = aTrack.GetMaterial();
222  G4VPhysicalVolume* aPVolume = aTrack.GetVolume();
223  G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
224 
225  // G4Material* bMaterial = aLVolume->GetMaterial();
226  G4PolarizationManager * polarizationManger = G4PolarizationManager::GetInstance();
227 
228  const G4bool VolumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
229  G4StokesVector ElectronPolarization = polarizationManger->GetVolumePolarization(aLVolume);
230 
231  if (VolumeIsPolarized) {
232 
233  if (verboseLevel>=2) {
234  G4cout << "G4PolarizedCompton::ComputeSaturationFactor: " << G4endl;
235  G4cout << " Mom " << GammaDirection0 << G4endl;
236  G4cout << " Polarization " << GammaPolarization << G4endl;
237  G4cout << " MaterialPol. " << ElectronPolarization << G4endl;
238  G4cout << " Phys. Volume " << aPVolume->GetName() << G4endl;
239  G4cout << " Log. Volume " << aLVolume->GetName() << G4endl;
240  G4cout << " Material " << aMaterial << G4endl;
241  }
242 
243  size_t midx = CurrentMaterialCutsCoupleIndex();
244  const G4PhysicsVector* aVector = nullptr;
245  if(midx < theAsymmetryTable->size()) {
246  aVector = (*theAsymmetryTable)(midx);
247  }
248  if (aVector) {
249  G4double asymmetry = aVector->Value(GammaEnergy);
250 
251  // we have to determine angle between particle motion
252  // and target polarisation here
253  // circ pol * Vec(ElectronPol)*Vec(PhotonMomentum)
254  // both vectors in global reference frame
255 
256  G4double pol = ElectronPolarization*GammaDirection0;
257  G4double polProduct = GammaPolarization.p3() * pol;
258  factor /= (1. + polProduct * asymmetry);
259  if (verboseLevel>=2) {
260  G4cout << " Asymmetry: " << asymmetry << G4endl;
261  G4cout << " PolProduct: " << polProduct << G4endl;
262  G4cout << " Factor: " << factor << G4endl;
263  }
264  } else {
265  G4ExceptionDescription ed;
266  ed << "Problem with asymmetry table: material index " << midx
267  << " is out of range or the table is not filled";
268  G4Exception("G4PolarizedComptonModel::ComputeSaturationFactor","em0048",
269  JustWarning, ed, "");
270  }
271  }
272  return factor;
273 }
274 
275 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
276 
277 void G4PolarizedCompton::BuildPhysicsTable(const G4ParticleDefinition& part)
278 {
279  // *** build (unpolarized) cross section tables (Lambda)
280  G4VEmProcess::BuildPhysicsTable(part);
281  if(buildAsymmetryTable && emModel) {
282  G4bool isMaster = true;
283  const G4PolarizedCompton* masterProcess =
284  static_cast<const G4PolarizedCompton*>(GetMasterProcess());
285  if(masterProcess && masterProcess != this) { isMaster = false; }
286  if(isMaster) { BuildAsymmetryTable(part); }
287  }
288 }
289 
290 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
291 
292 void G4PolarizedCompton::BuildAsymmetryTable(const G4ParticleDefinition& part)
293 {
294  // cleanup old, initialise new table
295  CleanTable();
296  theAsymmetryTable =
297  G4PhysicsTableHelper::PreparePhysicsTable(theAsymmetryTable);
298 
299  // Access to materials
300  const G4ProductionCutsTable* theCoupleTable=
301  G4ProductionCutsTable::GetProductionCutsTable();
302  size_t numOfCouples = theCoupleTable->GetTableSize();
303  if(!theAsymmetryTable) { return; }
304  G4int nbins = LambdaBinning();
305  G4double emin = MinKinEnergy();
306  G4double emax = MaxKinEnergy();
307  G4PhysicsLogVector* aVector = nullptr;
308  G4PhysicsLogVector* bVector = nullptr;
309 
310  for(size_t i=0; i<numOfCouples; ++i) {
311  if (theAsymmetryTable->GetFlag(i)) {
312 
313  // create physics vector and fill it
314  const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
315  // use same parameters as for lambda
316  if(!aVector) {
317  aVector = new G4PhysicsLogVector(emin, emax, nbins);
318  aVector->SetSpline(true);
319  bVector = aVector;
320  } else {
321  bVector = new G4PhysicsLogVector(*aVector);
322  }
323 
324  for (G4int j = 0; j <= nbins; ++j ) {
325  G4double energy = bVector->Energy(j);
326  G4double tasm=0.;
327  G4double asym = ComputeAsymmetry(energy, couple, part, 0., tasm);
328  bVector->PutValue(j,asym);
329  }
330  G4PhysicsTableHelper::SetPhysicsVector(theAsymmetryTable, i, bVector);
331  }
332  }
333 }
334 
335 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
336 
337 G4double G4PolarizedCompton::ComputeAsymmetry(G4double energy,
338  const G4MaterialCutsCouple* couple,
339  const G4ParticleDefinition& aParticle,
340  G4double cut,
341  G4double & tAsymmetry)
342 {
343  G4double lAsymmetry = 0.0;
344  tAsymmetry=0;
345 
346  //
347  // calculate polarized cross section
348  //
349  G4ThreeVector thePolarization=G4ThreeVector(0.,0.,1.);
350  emModel->SetTargetPolarization(thePolarization);
351  emModel->SetBeamPolarization(thePolarization);
352  G4double sigma2=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
353 
354  //
355  // calculate unpolarized cross section
356  //
357  thePolarization=G4ThreeVector();
358  emModel->SetTargetPolarization(thePolarization);
359  emModel->SetBeamPolarization(thePolarization);
360  G4double sigma0=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
361 
362  // determine assymmetries
363  if (sigma0 > 0.) {
364  lAsymmetry = sigma2/sigma0-1.;
365  }
366  return lAsymmetry;
367 }
368 
369 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....