Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ePolarizedIonisation.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: G4ePolarizedIonisation.cc 97175 2016-05-27 12:42:05Z gcosmo $
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4ePolarizedIonisation
34 //
35 // Author: A.Schaelicke on base of Vladimir Ivanchenko code
36 //
37 // Creation date: 10.11.2005
38 //
39 // Modifications:
40 //
41 // 10-11-05, include polarization description (A.Schaelicke)
42 // , create asymmetry table and determine interactionlength
43 // , update polarized differential cross section
44 //
45 // 20-08-06, modified interface (A.Schaelicke)
46 // 11-06-07, add PostStepGetPhysicalInteractionLength (A.Schalicke)
47 //
48 // Class Description:
49 //
50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
52 
54 #include "G4Electron.hh"
56 #include "G4UnitsTable.hh"
57 
59 #include "G4PhysicsTableHelper.hh"
60 #include "G4ProductionCutsTable.hh"
61 #include "G4PolarizationManager.hh"
62 #include "G4PolarizationHelper.hh"
63 #include "G4StokesVector.hh"
64 #include "G4EmParameters.hh"
65 
66 #include "G4SystemOfUnits.hh"
67 
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
69 
71  : G4VEnergyLossProcess(name),
72  theElectron(G4Electron::Electron()),
73  isElectron(true),
74  isInitialised(false),
75  theTargetPolarization(0.,0.,0.),
76  theAsymmetryTable(nullptr),
77  theTransverseAsymmetryTable(nullptr)
78 {
79  verboseLevel=0;
81  SetSecondaryParticle(theElectron);
82  flucModel = nullptr;
83  emModel = nullptr;
84 }
85 
86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
87 
89 {
90  CleanTables();
91 }
92 
93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
94 
95 void G4ePolarizedIonisation::CleanTables()
96 {
97  if(theAsymmetryTable) {
98  theAsymmetryTable->clearAndDestroy();
99  delete theAsymmetryTable;
100  theAsymmetryTable = nullptr;
101  }
102  if(theTransverseAsymmetryTable) {
103  theTransverseAsymmetryTable->clearAndDestroy();
104  delete theTransverseAsymmetryTable;
105  theTransverseAsymmetryTable = nullptr;
106  }
107 }
108 
109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
110 
111 G4double
113  const G4Material*, G4double cut)
114 {
115  G4double x = cut;
116  if(isElectron) { x += cut; }
117  return x;
118 }
119 
120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
121 
123 {
124  return (&p == G4Electron::Electron() || &p == G4Positron::Positron());
125 }
126 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
127 
129  const G4ParticleDefinition* part,
130  const G4ParticleDefinition* /*part2*/)
131 {
132  if(!isInitialised) {
133 
134  if(part == G4Positron::Positron()) { isElectron = false; }
135 
137  flucModel = FluctModel();
138 
139  emModel = new G4PolarizedMollerBhabhaModel();
140  SetEmModel(emModel, 1);
142  emModel->SetLowEnergyLimit(param->MinKinEnergy());
143  emModel->SetHighEnergyLimit(param->MaxKinEnergy());
144  AddEmModel(1, emModel, flucModel);
145 
146  isInitialised = true;
147  }
148 }
149 
150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
151 
153 {}
154 
155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
156 
158  G4double step,
159  G4ForceCondition* cond)
160 {
161  // *** get unploarised mean free path from lambda table ***
162  G4double mfp = G4VEnergyLossProcess::GetMeanFreePath(track, step, cond);
163  if(theAsymmetryTable && theTransverseAsymmetryTable && mfp < DBL_MAX) {
164  mfp *= ComputeSaturationFactor(track);
165  }
166  if (verboseLevel>=2) {
167  G4cout << "G4ePolarizedIonisation::MeanFreePath: "
168  << mfp / mm << " mm " << G4endl;
169  }
170  return mfp;
171 }
172 
173 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
174 
176  G4double step,
177  G4ForceCondition* cond)
178 {
179  // save previous values
182 
183  // *** get unpolarised mean free path from lambda table ***
184  // this changes theNumberOfInteractionLengthLeft and currentInteractionLength
186  G4double x0 = x;
187  G4double satFact = 1.;
188 
189  // *** add corrections on polarisation ***
190  if(theAsymmetryTable && theTransverseAsymmetryTable && x < DBL_MAX) {
191  satFact = ComputeSaturationFactor(track);
192  G4double curLength = currentInteractionLength*satFact;
193  G4double prvLength = iLength*satFact;
194  if(nLength > 0.0) {
196  std::max(nLength - step/prvLength, 0.0);
197  }
198  x = theNumberOfInteractionLengthLeft * curLength;
199  }
200  if (verboseLevel>=2) {
201  G4cout << "G4ePolarizedIonisation::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
212 G4ePolarizedIonisation::ComputeSaturationFactor(const G4Track& track)
213 {
214  G4Material* aMaterial = track.GetMaterial();
215  G4VPhysicalVolume* aPVolume = track.GetVolume();
216  G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
217 
219 
220  const G4bool volumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
221  G4StokesVector volPolarization = polarizationManger->GetVolumePolarization(aLVolume);
222 
223  G4double factor = 1.0;
224 
225  if (volumeIsPolarized && !volPolarization.IsZero()) {
226 
227  // *** get asymmetry, if target is polarized ***
228  const G4DynamicParticle* aDynamicPart = track.GetDynamicParticle();
229  const G4double energy = aDynamicPart->GetKineticEnergy();
230  const G4StokesVector polarization = track.GetPolarization();
231  const G4ParticleMomentum direction0 = aDynamicPart->GetMomentumDirection();
232 
233  if (verboseLevel>=2) {
234  G4cout << "G4ePolarizedIonisation::ComputeSaturationFactor: " << G4endl;
235  G4cout << " Energy(MeV) " << energy/MeV << G4endl;
236  G4cout << " Direction " << direction0 << G4endl;
237  G4cout << " Polarization " << polarization << G4endl;
238  G4cout << " MaterialPol. " << volPolarization << 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  const G4PhysicsVector* bVector = nullptr;
247  if(midx < theAsymmetryTable->size()) {
248  aVector = (*theAsymmetryTable)(midx);
249  }
250  if(midx < theTransverseAsymmetryTable->size()) {
251  bVector = (*theTransverseAsymmetryTable)(midx);
252  }
253  if(aVector && bVector) {
254  G4double lAsymmetry = aVector->Value(energy);
255  G4double tAsymmetry = bVector->Value(energy);
256  G4double polZZ = polarization.z()*(volPolarization*direction0);
257  G4double polXX = polarization.x()*
258  (volPolarization*G4PolarizationHelper::GetParticleFrameX(direction0));
259  G4double polYY = polarization.y()*
260  (volPolarization*G4PolarizationHelper::GetParticleFrameY(direction0));
261 
262  factor /= (1. + polZZ*lAsymmetry + (polXX + polYY)*tAsymmetry);
263 
264  if (verboseLevel>=2) {
265  G4cout << " Asymmetry: " << lAsymmetry << ", " << tAsymmetry << G4endl;
266  G4cout << " PolProduct: " << polXX << ", " << polYY << ", " << polZZ << G4endl;
267  G4cout << " Factor: " << factor << G4endl;
268  }
269  } else {
271  ed << "Problem with asymmetry tables: material index " << midx
272  << " is out of range or tables are not filled";
273  G4Exception("G4ePolarizedIonisation::ComputeSaturationFactor","em0048",
274  JustWarning, ed, "");
275  }
276  }
277  return factor;
278 }
279 
280 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
281 
283  const G4ParticleDefinition& part)
284 {
285  // *** build DEDX and (unpolarized) cross section tables
287  G4bool master = true;
288  const G4ePolarizedIonisation* masterProcess =
289  static_cast<const G4ePolarizedIonisation*>(GetMasterProcess());
290  if(masterProcess && masterProcess != this) { master = false; }
291  if(master) { BuildAsymmetryTables(part); }
292 }
293 
294 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
295 
296 void G4ePolarizedIonisation::BuildAsymmetryTables(
297  const G4ParticleDefinition& part)
298 {
299  // cleanup old, initialise new table
300  CleanTables();
301  theAsymmetryTable =
303  theTransverseAsymmetryTable =
304  G4PhysicsTableHelper::PreparePhysicsTable(theTransverseAsymmetryTable);
305 
306  const G4ProductionCutsTable* theCoupleTable=
308  size_t numOfCouples = theCoupleTable->GetTableSize();
309 
310  for (size_t j=0 ; j < numOfCouples; j++ ) {
311  // get cut value
312  const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
313 
314  G4double cut = (*theCoupleTable->GetEnergyCutsVector(1))[j];
315 
316  //create physics vectors then fill it (same parameters as lambda vector)
317  G4PhysicsVector * ptrVectorA = LambdaPhysicsVector(couple,cut);
318  G4PhysicsVector * ptrVectorB = LambdaPhysicsVector(couple,cut);
319  size_t bins = ptrVectorA->GetVectorLength();
320 
321  for (size_t i = 0 ; i < bins ; i++ ) {
322  G4double lowEdgeEnergy = ptrVectorA->Energy(i);
323  G4double tasm=0.;
324  G4double asym = ComputeAsymmetry(lowEdgeEnergy, couple, part, cut, tasm);
325  ptrVectorA->PutValue(i,asym);
326  ptrVectorB->PutValue(i,tasm);
327  }
328  theAsymmetryTable->insertAt( j , ptrVectorA ) ;
329  theTransverseAsymmetryTable->insertAt( j , ptrVectorB ) ;
330  }
331 }
332 
333 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
334 
335 G4double
336 G4ePolarizedIonisation::ComputeAsymmetry(G4double energy,
337  const G4MaterialCutsCouple* couple,
338  const G4ParticleDefinition& aParticle,
339  G4double cut,
340  G4double & tAsymmetry)
341 {
342  G4double lAsymmetry = 0.0;
343  tAsymmetry = 0.0;
344  if (isElectron) { lAsymmetry = tAsymmetry = -1.0; }
345 
346  // calculate polarized cross section
347  theTargetPolarization=G4ThreeVector(0.,0.,1.);
348  emModel->SetTargetPolarization(theTargetPolarization);
349  emModel->SetBeamPolarization(theTargetPolarization);
350  G4double sigma2=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
351 
352  // calculate transversely polarized cross section
353  theTargetPolarization=G4ThreeVector(1.,0.,0.);
354  emModel->SetTargetPolarization(theTargetPolarization);
355  emModel->SetBeamPolarization(theTargetPolarization);
356  G4double sigma3=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
357 
358  // calculate unpolarized cross section
359  theTargetPolarization=G4ThreeVector();
360  emModel->SetTargetPolarization(theTargetPolarization);
361  emModel->SetBeamPolarization(theTargetPolarization);
362  G4double sigma0 = emModel->CrossSection(couple,&aParticle,energy,cut,energy);
363  // determine assymmetries
364  if (sigma0 > 0.) {
365  lAsymmetry=sigma2/sigma0 - 1.;
366  tAsymmetry=sigma3/sigma0 - 1.;
367  }
368  if (std::fabs(lAsymmetry)>1.) {
369  G4cout<<"G4ePolarizedIonisation::ComputeAsymmetry WARNING: E(MeV)= "
370  << energy << " lAsymmetry= "<<lAsymmetry
371  <<" ("<<std::fabs(lAsymmetry)-1.<<")\n";
372  }
373  if (std::fabs(tAsymmetry)>1.) {
374  G4cout<<" energy="<<energy<<"\n";
375  G4cout<<"G4ePolarizedIonisation::ComputeAsymmetry WARNING: E(MeV)= "
376  << energy << " tAsymmetry= "<<tAsymmetry
377  <<" ("<<std::fabs(tAsymmetry)-1.<<")\n";
378  }
379  return lAsymmetry;
380 }
381 
382 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
383 
const G4VProcess * GetMasterProcess() const
Definition: G4VProcess.hh:538
const XML_Char * name
Definition: expat.h:151
const G4ThreeVector & GetPolarization() const
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
G4double MaxKinEnergy() const
static constexpr double mm
Definition: G4SIunits.hh:115
G4bool isElectron(G4int ityp)
G4int verboseLevel
Definition: G4VProcess.hh:368
G4ePolarizedIonisation(const G4String &name="pol-eIoni")
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
virtual void BuildPhysicsTable(const G4ParticleDefinition &) override
double x() const
const G4DynamicParticle * GetDynamicParticle() const
virtual G4bool IsApplicable(const G4ParticleDefinition &p) override
const char * p
Definition: xmltok.h:285
static G4PolarizationManager * GetInstance()
virtual G4double GetMeanFreePath(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition) override
virtual G4double GetMeanFreePath(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition) override
G4double theNumberOfInteractionLengthLeft
Definition: G4VProcess.hh:293
tuple x
Definition: test.py:50
size_t GetVectorLength() const
void SetFluctModel(G4VEmFluctuationModel *)
G4PhysicsVector * LambdaPhysicsVector(const G4MaterialCutsCouple *, G4double cut)
static G4PhysicsTable * PreparePhysicsTable(G4PhysicsTable *physTable)
double z() const
virtual G4double MinPrimaryEnergy(const G4ParticleDefinition *, const G4Material *, G4double cut) override
G4bool IsZero() const
virtual G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition) override
virtual void InitialiseEnergyLossProcess(const G4ParticleDefinition *, const G4ParticleDefinition *) override
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:732
void AddEmModel(G4int, G4VEmModel *, G4VEmFluctuationModel *fluc=0, const G4Region *region=nullptr)
G4VEmFluctuationModel * FluctModel()
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:503
G4double currentInteractionLength
Definition: G4VProcess.hh:297
void PutValue(size_t index, G4double theValue)
virtual void PrintInfo() override
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
G4double MinKinEnergy() const
G4double Energy(size_t index) const
G4double Value(G4double theEnergy, size_t &lastidx) const
void SetSecondaryParticle(const G4ParticleDefinition *p)
G4Material * GetMaterial() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4ProductionCutsTable * GetProductionCutsTable()
static G4Positron * Positron()
Definition: G4Positron.cc:94
virtual void BuildPhysicsTable(const G4ParticleDefinition &) override
G4LogicalVolume * GetLogicalVolume() const
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double energy(const ThreeVector &p, const G4double m)
static G4EmParameters * Instance()
double y() const
bool IsPolarized(G4LogicalVolume *lVol) const
void SetTargetPolarization(const G4ThreeVector &pTarget)
void SetEmModel(G4VEmModel *, G4int index=1)
G4VPhysicalVolume * GetVolume() const
static G4ThreeVector GetParticleFrameY(const G4ThreeVector &)
void insertAt(size_t, G4PhysicsVector *)
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
double G4double
Definition: G4Types.hh:76
void SetBeamPolarization(const G4ThreeVector &pBeam)
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:739
G4ForceCondition
#define DBL_MAX
Definition: templates.hh:83
static G4ThreeVector GetParticleFrameX(const G4ThreeVector &)
virtual G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition) override
const G4ThreeVector & GetVolumePolarization(G4LogicalVolume *lVol) const
void clearAndDestroy()
size_t CurrentMaterialCutsCoupleIndex() const