Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4PolarizedCompton.cc
Go to the documentation of this file.
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 // $Id: G4PolarizedCompton.cc 97175 2016-05-27 12:42:05Z gcosmo $
28 //
29 //
30 // File name: G4PolarizedCompton
31 //
32 // Author: Andreas Schaelicke
33 // based on code by Michel Maire / Vladimir IVANTCHENKO
34 // Class description
35 //
36 // modified version respecting media and beam polarization
37 // using the stokes formalism
38 //
39 // Creation date: 01.05.2005
40 //
41 // Modifications:
42 //
43 // 01-01-05, include polarization description (A.Stahl)
44 // 01-01-05, create asymmetry table and determine interactionlength (A.Stahl)
45 // 01-05-05, update handling of media polarization (A.Schalicke)
46 // 01-05-05, update polarized differential cross section (A.Schalicke)
47 // 20-05-05, added polarization transfer (A.Schalicke)
48 // 10-06-05, transformation between different reference frames (A.Schalicke)
49 // 17-10-05, correct reference frame dependence in GetMeanFreePath (A.Schalicke)
50 // 26-07-06, cross section recalculated (P.Starovoitov)
51 // 09-08-06, make it work under current geant4 release (A.Schalicke)
52 // 11-06-07, add PostStepGetPhysicalInteractionLength (A.Schalicke)
53 // -----------------------------------------------------------------------------
54 
55 
56 #include "G4PolarizedCompton.hh"
57 #include "G4SystemOfUnits.hh"
58 #include "G4Electron.hh"
59 
60 #include "G4StokesVector.hh"
61 #include "G4PolarizationManager.hh"
63 #include "G4ProductionCutsTable.hh"
64 #include "G4PhysicsTableHelper.hh"
65 #include "G4KleinNishinaCompton.hh"
67 #include "G4EmParameters.hh"
68 
69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
70 
71 G4PhysicsTable* G4PolarizedCompton::theAsymmetryTable = nullptr;
72 
74  G4ProcessType type):
75  G4VEmProcess (processName, type),
76  buildAsymmetryTable(true),
77  useAsymmetryTable(true),
78  isInitialised(false),
79  mType(10),
80  targetPolarization(0.0,0.0,0.0)
81 {
83  SetBuildTableFlag(true);
87  SetSplineFlag(true);
88  emModel = nullptr;
89 }
90 
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
92 
94 {
95  CleanTable();
96 }
97 
98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
99 
100 void G4PolarizedCompton::CleanTable()
101 {
102  if( theAsymmetryTable) {
103  theAsymmetryTable->clearAndDestroy();
104  delete theAsymmetryTable;
105  theAsymmetryTable = nullptr;
106  }
107 }
108 
109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
110 
112 {
113  return (&p == G4Gamma::Gamma());
114 }
115 
116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
117 
119 {
120  if(!isInitialised) {
121  isInitialised = true;
122  if(0 == mType) {
123  if(!EmModel(1)) { SetEmModel(new G4KleinNishinaCompton(), 1); }
124  } else {
125  emModel = new G4PolarizedComptonModel();
126  SetEmModel(emModel, 1);
127  }
129  EmModel(1)->SetLowEnergyLimit(param->MinKinEnergy());
131  AddEmModel(1, EmModel(1));
132  }
133 }
134 
135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
136 
138 {
139  G4cout << " Total cross sections has a good parametrisation"
140  << " from 10 KeV to (100/Z) GeV"
141  << "\n Sampling according " << EmModel(1)->GetName() << " model"
142  << G4endl;
143 }
144 
145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
146 
148 {
149  if(ss == "Klein-Nishina") { mType = 0; }
150  if(ss == "Polarized-Compton") { mType = 10; }
151 }
152 
153 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
154 
156  G4double previousStepSize,
158 {
159  // *** get unploarised mean free path from lambda table ***
160  G4double mfp = G4VEmProcess::GetMeanFreePath(aTrack, previousStepSize, condition);
161 
162  if (theAsymmetryTable && useAsymmetryTable && mfp < DBL_MAX) {
163  mfp *= ComputeSaturationFactor(aTrack);
164  }
165  if (verboseLevel>=2) {
166  G4cout << "G4PolarizedCompton::MeanFreePath: " << mfp / mm << " mm " << G4endl;
167  }
168  return mfp;
169 }
170 
171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
172 
174  const G4Track& aTrack,
175  G4double previousStepSize,
177 {
178  // save previous values
181 
182  // *** compute unpolarized step limit ***
183  // this changes theNumberOfInteractionLengthLeft and currentInteractionLength
185  previousStepSize,
186  condition);
187  G4double x0 = x;
188  G4double satFact = 1.0;
189 
190  // *** add corrections on polarisation ***
191  if (theAsymmetryTable && useAsymmetryTable && x < DBL_MAX) {
192  satFact = ComputeSaturationFactor(aTrack);
193  G4double curLength = currentInteractionLength*satFact;
194  G4double prvLength = iLength*satFact;
195  if(nLength > 0.0) {
197  std::max(nLength - previousStepSize/prvLength, 0.0);
198  }
199  x = theNumberOfInteractionLengthLeft * curLength;
200  }
201  if (verboseLevel>=2) {
202  G4cout << "G4PolarizedCompton::PostStepGPIL: "
203  << std::setprecision(8) << x/mm << " mm;" << G4endl
204  << " unpolarized value: "
205  << std::setprecision(8) << x0/mm << " mm." << G4endl;
206  }
207  return x;
208 }
209 
210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
211 
212 G4double G4PolarizedCompton::ComputeSaturationFactor(const G4Track& aTrack)
213 {
214  G4double factor = 1.0;
215 
216  // *** get asymmetry, if target is polarized ***
217  const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
218  const G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
219  const G4StokesVector GammaPolarization = aTrack.GetPolarization();
220  const G4ParticleMomentum GammaDirection0 = aDynamicGamma->GetMomentumDirection();
221 
222  G4Material* aMaterial = aTrack.GetMaterial();
223  G4VPhysicalVolume* aPVolume = aTrack.GetVolume();
224  G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
225 
226  // G4Material* bMaterial = aLVolume->GetMaterial();
228 
229  const G4bool VolumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
230  G4StokesVector ElectronPolarization = polarizationManger->GetVolumePolarization(aLVolume);
231 
232  if (VolumeIsPolarized) {
233 
234  if (verboseLevel>=2) {
235  G4cout << "G4PolarizedCompton::ComputeSaturationFactor: " << G4endl;
236  G4cout << " Mom " << GammaDirection0 << G4endl;
237  G4cout << " Polarization " << GammaPolarization << G4endl;
238  G4cout << " MaterialPol. " << ElectronPolarization << G4endl;
239  G4cout << " Phys. Volume " << aPVolume->GetName() << G4endl;
240  G4cout << " Log. Volume " << aLVolume->GetName() << G4endl;
241  G4cout << " Material " << aMaterial << G4endl;
242  }
243 
244  size_t midx = CurrentMaterialCutsCoupleIndex();
245  const G4PhysicsVector* aVector = nullptr;
246  if(midx < theAsymmetryTable->size()) {
247  aVector = (*theAsymmetryTable)(midx);
248  }
249  if (aVector) {
250  G4double asymmetry = aVector->Value(GammaEnergy);
251 
252  // we have to determine angle between particle motion
253  // and target polarisation here
254  // circ pol * Vec(ElectronPol)*Vec(PhotonMomentum)
255  // both vectors in global reference frame
256 
257  G4double pol = ElectronPolarization*GammaDirection0;
258  G4double polProduct = GammaPolarization.p3() * pol;
259  factor /= (1. + polProduct * asymmetry);
260  if (verboseLevel>=2) {
261  G4cout << " Asymmetry: " << asymmetry << G4endl;
262  G4cout << " PolProduct: " << polProduct << G4endl;
263  G4cout << " Factor: " << factor << G4endl;
264  }
265  } else {
267  ed << "Problem with asymmetry table: material index " << midx
268  << " is out of range or the table is not filled";
269  G4Exception("G4PolarizedComptonModel::ComputeSaturationFactor","em0048",
270  JustWarning, ed, "");
271  }
272  }
273  return factor;
274 }
275 
276 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
277 
279 {
280  // *** build (unpolarized) cross section tables (Lambda)
282  if(buildAsymmetryTable && emModel) {
283  G4bool isMaster = true;
284  const G4PolarizedCompton* masterProcess =
285  static_cast<const G4PolarizedCompton*>(GetMasterProcess());
286  if(masterProcess && masterProcess != this) { isMaster = false; }
287  if(isMaster) { BuildAsymmetryTable(part); }
288  }
289 }
290 
291 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
292 
293 void G4PolarizedCompton::BuildAsymmetryTable(const G4ParticleDefinition& part)
294 {
295  // cleanup old, initialise new table
296  CleanTable();
297  theAsymmetryTable =
299 
300  // Access to materials
301  const G4ProductionCutsTable* theCoupleTable=
303  size_t numOfCouples = theCoupleTable->GetTableSize();
304  if(!theAsymmetryTable) { return; }
305  G4int nbins = LambdaBinning();
306  G4double emin = MinKinEnergy();
308  G4PhysicsLogVector* aVector = nullptr;
309  G4PhysicsLogVector* bVector = nullptr;
310 
311  for(size_t i=0; i<numOfCouples; ++i) {
312  if (theAsymmetryTable->GetFlag(i)) {
313 
314  // create physics vector and fill it
315  const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
316  // use same parameters as for lambda
317  if(!aVector) {
318  aVector = new G4PhysicsLogVector(emin, emax, nbins);
319  aVector->SetSpline(true);
320  bVector = aVector;
321  } else {
322  bVector = new G4PhysicsLogVector(*aVector);
323  }
324 
325  for (G4int j = 0; j <= nbins; ++j ) {
326  G4double energy = bVector->Energy(j);
327  G4double tasm=0.;
328  G4double asym = ComputeAsymmetry(energy, couple, part, 0., tasm);
329  bVector->PutValue(j,asym);
330  }
331  G4PhysicsTableHelper::SetPhysicsVector(theAsymmetryTable, i, bVector);
332  }
333  }
334 }
335 
336 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
337 
338 G4double G4PolarizedCompton::ComputeAsymmetry(G4double energy,
339  const G4MaterialCutsCouple* couple,
340  const G4ParticleDefinition& aParticle,
341  G4double cut,
342  G4double & tAsymmetry)
343 {
344  G4double lAsymmetry = 0.0;
345  tAsymmetry=0;
346 
347  //
348  // calculate polarized cross section
349  //
350  G4ThreeVector thePolarization=G4ThreeVector(0.,0.,1.);
351  emModel->SetTargetPolarization(thePolarization);
352  emModel->SetBeamPolarization(thePolarization);
353  G4double sigma2=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
354 
355  //
356  // calculate unpolarized cross section
357  //
358  thePolarization=G4ThreeVector();
359  emModel->SetTargetPolarization(thePolarization);
360  emModel->SetBeamPolarization(thePolarization);
361  G4double sigma0=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
362 
363  // determine assymmetries
364  if (sigma0 > 0.) {
365  lAsymmetry = sigma2/sigma0-1.;
366  }
367  return lAsymmetry;
368 }
369 
370 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
G4double condition(const G4ErrorSymMatrix &m)
const G4VProcess * GetMasterProcess() const
Definition: G4VProcess.hh:538
const G4ThreeVector & GetPolarization() const
static void SetPhysicsVector(G4PhysicsTable *physTable, size_t idx, G4PhysicsVector *vec)
virtual void BuildPhysicsTable(const G4ParticleDefinition &) override
G4double MaxKinEnergy() const
static constexpr double mm
Definition: G4SIunits.hh:115
virtual void BuildPhysicsTable(const G4ParticleDefinition &) override
G4int verboseLevel
Definition: G4VProcess.hh:368
virtual G4bool IsApplicable(const G4ParticleDefinition &) override
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
const G4DynamicParticle * GetDynamicParticle() const
void SetBuildTableFlag(G4bool val)
G4VEmModel * EmModel(G4int index=1) const
const char * p
Definition: xmltok.h:285
static G4PolarizationManager * GetInstance()
virtual G4double GetMeanFreePath(const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition) override
void SetSplineFlag(G4bool val)
G4double MaxKinEnergy() const
G4double theNumberOfInteractionLengthLeft
Definition: G4VProcess.hh:293
virtual void PrintInfo() override
void SetStartFromNullFlag(G4bool val)
virtual G4double GetMeanFreePath(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition) override
int G4int
Definition: G4Types.hh:78
static G4PhysicsTable * PreparePhysicsTable(G4PhysicsTable *physTable)
virtual void InitialiseProcess(const G4ParticleDefinition *) override
G4int LambdaBinning() const
void SetSpline(G4bool)
G4double p3() const
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:724
void SetMinKinEnergyPrim(G4double e)
void SetEmModel(G4VEmModel *, G4int index=1)
size_t CurrentMaterialCutsCoupleIndex() const
G4GLOB_DLL std::ostream G4cout
const G4String & GetName() const
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
G4double CrossSection(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:500
G4double currentInteractionLength
Definition: G4VProcess.hh:297
void PutValue(size_t index, G4double theValue)
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
virtual G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition) override
G4double MinKinEnergy() const
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4double Energy(size_t index) const
G4double Value(G4double theEnergy, size_t &lastidx) const
G4Material * GetMaterial() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static const G4double emax
void AddEmModel(G4int, G4VEmModel *, const G4Region *region=nullptr)
static G4ProductionCutsTable * GetProductionCutsTable()
G4PolarizedCompton(const G4String &processName="pol-compt", G4ProcessType type=fElectromagnetic)
G4LogicalVolume * GetLogicalVolume() const
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
void SetBeamPolarization(const G4ThreeVector &pBeam)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double energy(const ThreeVector &p, const G4double m)
void SetSecondaryParticle(const G4ParticleDefinition *p)
static G4EmParameters * Instance()
bool IsPolarized(G4LogicalVolume *lVol) const
G4VPhysicalVolume * GetVolume() const
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
const G4String & GetName() const
Definition: G4VEmModel.hh:794
G4bool GetFlag(size_t i) const
const G4String & GetName() const
double G4double
Definition: G4Types.hh:76
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:731
G4ForceCondition
G4double MinKinEnergy() const
#define DBL_MAX
Definition: templates.hh:83
const G4ThreeVector & GetVolumePolarization(G4LogicalVolume *lVol) const
void SetTargetPolarization(const G4ThreeVector &pTarget)
void clearAndDestroy()
void SetModel(const G4String &name)
virtual G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition) override
G4ProcessType