Geant4  10.02.p01
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 93113 2015-10-07 07:49:04Z 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 
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 
101 {
102  if( theAsymmetryTable) {
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 {
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 value
180 
181  // *** compute uppolarized step limit ***
183  previousStepSize,
184  condition);
185 
186  // *** add corrections on polarisation ***
189  if(nLength > 0.0) {
191  std::max(nLength - previousStepSize/curLength, 0.0);
192  }
193  x = theNumberOfInteractionLengthLeft * curLength;
194  }
195  if (verboseLevel>=2) {
196  G4cout << "G4PolarizedCompton::PostStepGetPhysicalInteractionLength: "
197  << x/mm << " mm " << G4endl;
198  }
199  return x;
200 }
201 
202 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
203 
205 {
206  G4double factor = 1.0;
207 
208  // *** get asymmetry, if target is polarized ***
209  const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
210  const G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
211  const G4StokesVector GammaPolarization = aTrack.GetPolarization();
212  const G4ParticleMomentum GammaDirection0 = aDynamicGamma->GetMomentumDirection();
213 
214  G4Material* aMaterial = aTrack.GetMaterial();
215  G4VPhysicalVolume* aPVolume = aTrack.GetVolume();
216  G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
217 
218  // G4Material* bMaterial = aLVolume->GetMaterial();
220 
221  const G4bool VolumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
222  G4StokesVector ElectronPolarization = polarizationManger->GetVolumePolarization(aLVolume);
223 
224  if (VolumeIsPolarized) {
225 
226  if (verboseLevel>=2) {
227  G4cout << "G4PolarizedCompton::ComputeSaturationFactor: " << G4endl;
228  G4cout << " Mom " << GammaDirection0 << G4endl;
229  G4cout << " Polarization " << GammaPolarization << G4endl;
230  G4cout << " MaterialPol. " << ElectronPolarization << G4endl;
231  G4cout << " Phys. Volume " << aPVolume->GetName() << G4endl;
232  G4cout << " Log. Volume " << aLVolume->GetName() << G4endl;
233  G4cout << " Material " << aMaterial << G4endl;
234  }
235 
236  size_t midx = CurrentMaterialCutsCoupleIndex();
237  const G4PhysicsVector* aVector = nullptr;
238  if(midx < theAsymmetryTable->size()) {
239  aVector = (*theAsymmetryTable)(midx);
240  }
241  if (aVector) {
242  G4double asymmetry = aVector->Value(GammaEnergy);
243 
244  // we have to determine angle between particle motion
245  // and target polarisation here
246  // circ pol * Vec(ElectronPol)*Vec(PhotonMomentum)
247  // both vectors in global reference frame
248 
249  G4double pol = ElectronPolarization*GammaDirection0;
250  G4double polProduct = GammaPolarization.p3() * pol;
251  factor /= (1. + polProduct * asymmetry);
252  if (verboseLevel>=2) {
253  G4cout << " Asymmetry: " << asymmetry << G4endl;
254  G4cout << " PolProduct: " << polProduct << G4endl;
255  G4cout << " Factor: " << factor << G4endl;
256  }
257  } else {
259  ed << "Problem with asymmetry table: material index " << midx
260  << " is out of range or the table is not filled";
261  G4Exception("G4PolarizedComptonModel::ComputeSaturationFactor","em0048",
262  JustWarning, ed, "");
263  }
264  }
265  return factor;
266 }
267 
268 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
269 
271 {
272  // *** build (unpolarized) cross section tables (Lambda)
274  if(buildAsymmetryTable && emModel) {
275  G4bool isMaster = true;
276  const G4PolarizedCompton* masterProcess =
277  static_cast<const G4PolarizedCompton*>(GetMasterProcess());
278  if(masterProcess && masterProcess != this) { isMaster = false; }
279  if(isMaster) { BuildAsymmetryTable(part); }
280  }
281 }
282 
283 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
284 
286 {
287  // cleanup old, initialise new table
288  CleanTable();
291 
292  // Access to materials
293  const G4ProductionCutsTable* theCoupleTable=
295  size_t numOfCouples = theCoupleTable->GetTableSize();
296  if(!theAsymmetryTable) { return; }
297  G4int nbins = LambdaBinning();
298  G4double emin = MinKinEnergy();
300  G4PhysicsLogVector* aVector = 0;
301  G4PhysicsLogVector* bVector = 0;
302 
303  for(size_t i=0; i<numOfCouples; ++i) {
304  if (theAsymmetryTable->GetFlag(i)) {
305 
306  // create physics vector and fill it
307  const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
308  // use same parameters as for lambda
309  if(!aVector) {
310  aVector = new G4PhysicsLogVector(emin, emax, nbins);
311  aVector->SetSpline(true);
312  bVector = aVector;
313  } else {
314  bVector = new G4PhysicsLogVector(*aVector);
315  }
316 
317  for (G4int j = 0; j <= nbins; ++j ) {
318  G4double energy = bVector->Energy(j);
319  G4double tasm=0.;
320  G4double asym = ComputeAsymmetry(energy, couple, part, 0., tasm);
321  bVector->PutValue(j,asym);
322  }
324  }
325  }
326 }
327 
328 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
329 
331  const G4MaterialCutsCouple* couple,
332  const G4ParticleDefinition& aParticle,
333  G4double cut,
334  G4double & tAsymmetry)
335 {
336  G4double lAsymmetry = 0.0;
337  tAsymmetry=0;
338 
339  //
340  // calculate polarized cross section
341  //
342  G4ThreeVector thePolarization=G4ThreeVector(0.,0.,1.);
343  emModel->SetTargetPolarization(thePolarization);
344  emModel->SetBeamPolarization(thePolarization);
345  G4double sigma2=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
346 
347  //
348  // calculate unpolarized cross section
349  //
350  thePolarization=G4ThreeVector();
351  emModel->SetTargetPolarization(thePolarization);
352  emModel->SetBeamPolarization(thePolarization);
353  G4double sigma0=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
354 
355  // determine assymmetries
356  if (sigma0 > 0.) {
357  lAsymmetry = sigma2/sigma0-1.;
358  }
359  return lAsymmetry;
360 }
361 
362 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
G4double condition(const G4ErrorSymMatrix &m)
const G4VProcess * GetMasterProcess() const
Definition: G4VProcess.hh:538
const G4ThreeVector & GetPolarization() const
static const double MeV
Definition: G4SIunits.hh:211
static void SetPhysicsVector(G4PhysicsTable *physTable, size_t idx, G4PhysicsVector *vec)
G4double MaxKinEnergy() const
G4int verboseLevel
Definition: G4VProcess.hh:368
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
const G4DynamicParticle * GetDynamicParticle() const
void SetBuildTableFlag(G4bool val)
void BuildAsymmetryTable(const G4ParticleDefinition &part)
G4VEmModel * EmModel(G4int index=1) const
static G4PolarizationManager * GetInstance()
void SetSplineFlag(G4bool val)
G4double MaxKinEnergy() const
G4double theNumberOfInteractionLengthLeft
Definition: G4VProcess.hh:293
G4PolarizedComptonModel * emModel
virtual void BuildPhysicsTable(const G4ParticleDefinition &)
virtual G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
virtual G4bool IsApplicable(const G4ParticleDefinition &)
void SetStartFromNullFlag(G4bool val)
int G4int
Definition: G4Types.hh:78
static G4PhysicsTable * PreparePhysicsTable(G4PhysicsTable *physTable)
G4int LambdaBinning() const
virtual void BuildPhysicsTable(const G4ParticleDefinition &)
void SetSpline(G4bool)
G4double p3() const
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:725
void SetMinKinEnergyPrim(G4double e)
G4double ComputeSaturationFactor(const G4Track &aTrack)
virtual G4double GetMeanFreePath(const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition)
void SetEmModel(G4VEmModel *, G4int index=1)
size_t CurrentMaterialCutsCoupleIndex() const
G4GLOB_DLL std::ostream G4cout
const G4String & GetName() const
virtual G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
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:501
G4double currentInteractionLength
Definition: G4VProcess.hh:297
void PutValue(size_t index, G4double theValue)
virtual G4double GetMeanFreePath(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
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
static G4ProductionCutsTable * GetProductionCutsTable()
G4PolarizedCompton(const G4String &processName="pol-compt", G4ProcessType type=fElectromagnetic)
G4double ComputeAsymmetry(G4double energy, const G4MaterialCutsCouple *couple, const G4ParticleDefinition &particle, G4double cut, G4double &tAsymmetry)
static const G4double factor
G4LogicalVolume * GetLogicalVolume() const
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
void SetBeamPolarization(const G4ThreeVector &pBeam)
void AddEmModel(G4int, G4VEmModel *, const G4Region *region=0)
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)
const G4double x[NPOINTSGL]
static G4EmParameters * Instance()
bool IsPolarized(G4LogicalVolume *lVol) const
G4VPhysicalVolume * GetVolume() const
virtual void InitialiseProcess(const G4ParticleDefinition *)
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
const G4String & GetName() const
Definition: G4VEmModel.hh:795
G4bool GetFlag(size_t i) const
const G4String & GetName() const
double G4double
Definition: G4Types.hh:76
static G4PhysicsTable * theAsymmetryTable
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:732
G4ForceCondition
G4ThreeVector G4ParticleMomentum
G4double MinKinEnergy() const
#define DBL_MAX
Definition: templates.hh:83
static const double mm
Definition: G4SIunits.hh:114
const G4ThreeVector & GetVolumePolarization(G4LogicalVolume *lVol) const
void SetTargetPolarization(const G4ThreeVector &pTarget)
virtual void PrintInfo()
void clearAndDestroy()
void SetModel(const G4String &name)
G4ProcessType