Geant4  10.01.p02
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 85018 2014-10-23 09:51:37Z 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 "G4BohrFluctuations.hh"
57 #include "G4UnitsTable.hh"
58 
60 #include "G4ProductionCutsTable.hh"
61 #include "G4PolarizationManager.hh"
62 #include "G4PolarizationHelper.hh"
63 #include "G4StokesVector.hh"
64 #include "G4EmParameters.hh"
65 
66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
67 
69  : G4VEnergyLossProcess(name),
70  theElectron(G4Electron::Electron()),
71  isElectron(true),
72  isInitialised(false),
73  theAsymmetryTable(NULL),
74  theTransverseAsymmetryTable(NULL)
75 {
76  verboseLevel=0;
79  flucModel = 0;
80  emModel = 0;
81 }
82 
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
84 
86 {
87  delete theAsymmetryTable;
89 }
90 
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92 
94  const G4ParticleDefinition* part,
95  const G4ParticleDefinition* /*part2*/)
96 {
97  if(!isInitialised) {
98 
99  if(part == G4Positron::Positron()) { isElectron = false; }
100 
102  flucModel = FluctModel();
103 
105  SetEmModel(emModel, 1);
110 
111  isInitialised = true;
112  }
113 }
114 
115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
116 
118 {}
119 
120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
121 
123  G4double step,
124  G4ForceCondition* cond)
125 {
126  // *** get unploarised mean free path from lambda table ***
127  G4double mfp = G4VEnergyLossProcess::GetMeanFreePath(track, step, cond);
128 
129 
130  // *** get asymmetry, if target is polarized ***
131  G4VPhysicalVolume* aPVolume = track.GetVolume();
132  G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
133 
135  G4bool volumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
136  const G4StokesVector ePolarization = track.GetPolarization();
137 
138  if (mfp != DBL_MAX && volumeIsPolarized && !ePolarization.IsZero()) {
139  const G4DynamicParticle* aDynamicElectron = track.GetDynamicParticle();
140  G4double eEnergy = aDynamicElectron->GetKineticEnergy();
141  const G4ParticleMomentum eDirection0 = aDynamicElectron->GetMomentumDirection();
142 
143  G4StokesVector volumePolarization = polarizationManger->GetVolumePolarization(aLVolume);
144 
145  G4bool isOutRange;
146  size_t idx = CurrentMaterialCutsCoupleIndex();
147  G4double lAsymmetry = (*theAsymmetryTable)(idx)->
148  GetValue(eEnergy, isOutRange);
149  G4double tAsymmetry = (*theTransverseAsymmetryTable)(idx)->
150  GetValue(eEnergy, isOutRange);
151 
152  // calculate longitudinal spin component
153  G4double polZZ = ePolarization.z()*
154  volumePolarization*eDirection0;
155  // calculate transvers spin components
156  G4double polXX = ePolarization.x()*
157  volumePolarization*G4PolarizationHelper::GetParticleFrameX(eDirection0);
158  G4double polYY = ePolarization.y()*
159  volumePolarization*G4PolarizationHelper::GetParticleFrameY(eDirection0);
160 
161 
162  G4double impact = 1. + polZZ*lAsymmetry + (polXX + polYY)*tAsymmetry;
163  // determine polarization dependent mean free path
164  mfp /= impact;
165  if (mfp <=0.) {
166  G4cout <<"PV impact ( "<<polXX<<" , "<<polYY<<" , "<<polZZ<<" )"<<G4endl;
167  G4cout << " impact on MFP is "<< impact <<G4endl;
168  G4cout<<" lAsymmetry= "<<lAsymmetry<<" ("<<std::fabs(lAsymmetry)-1.<<")\n";
169  G4cout<<" tAsymmetry= "<<tAsymmetry<<" ("<<std::fabs(tAsymmetry)-1.<<")\n";
170  }
171  }
172 
173  return mfp;
174 }
175 
177  G4double step,
178  G4ForceCondition* cond)
179 {
180  // *** get unploarised mean free path from lambda table ***
182 
183 
184  // *** get asymmetry, if target is polarized ***
185  G4VPhysicalVolume* aPVolume = track.GetVolume();
186  G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
187 
189  G4bool volumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
190  const G4StokesVector ePolarization = track.GetPolarization();
191 
192  if (mfp != DBL_MAX && volumeIsPolarized && !ePolarization.IsZero()) {
193  const G4DynamicParticle* aDynamicElectron = track.GetDynamicParticle();
194  G4double eEnergy = aDynamicElectron->GetKineticEnergy();
195  const G4ParticleMomentum eDirection0 = aDynamicElectron->GetMomentumDirection();
196 
197  G4StokesVector volumePolarization = polarizationManger->GetVolumePolarization(aLVolume);
198 
199  size_t idx = CurrentMaterialCutsCoupleIndex();
200  G4double lAsymmetry = (*theAsymmetryTable)(idx)->Value(eEnergy);
201  G4double tAsymmetry = (*theTransverseAsymmetryTable)(idx)->Value(eEnergy);
202 
203  // calculate longitudinal spin component
204  G4double polZZ = ePolarization.z()*
205  volumePolarization*eDirection0;
206  // calculate transvers spin components
207  G4double polXX = ePolarization.x()*
208  volumePolarization*G4PolarizationHelper::GetParticleFrameX(eDirection0);
209  G4double polYY = ePolarization.y()*
210  volumePolarization*G4PolarizationHelper::GetParticleFrameY(eDirection0);
211 
212 
213  G4double impact = 1. + polZZ*lAsymmetry + (polXX + polYY)*tAsymmetry;
214  // determine polarization dependent mean free path
215  mfp /= impact;
216  if (mfp <=0.) {
217  G4cout <<"PV impact ( "<<polXX<<" , "<<polYY<<" , "<<polZZ<<" )"<<G4endl;
218  G4cout << " impact on MFP is "<< impact <<G4endl;
219  G4cout<<" lAsymmetry= "<<lAsymmetry<<" ("<<std::fabs(lAsymmetry)-1.<<")\n";
220  G4cout<<" tAsymmetry= "<<tAsymmetry<<" ("<<std::fabs(tAsymmetry)-1.<<")\n";
221  }
222  }
223 
224  return mfp;
225 }
226 
227 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
229 {
230  // *** build DEDX and (unpolarized) cross section tables
232  // G4PhysicsTable* pt =
233  // BuildDEDXTable();
234 
235 
236  // *** build asymmetry-table
237  if (theAsymmetryTable) {
241 
242  const G4ProductionCutsTable* theCoupleTable=
244  size_t numOfCouples = theCoupleTable->GetTableSize();
245 
246  theAsymmetryTable = new G4PhysicsTable(numOfCouples);
247  theTransverseAsymmetryTable = new G4PhysicsTable(numOfCouples);
248 
249  for (size_t j=0 ; j < numOfCouples; j++ ) {
250  // get cut value
251  const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
252 
253  G4double cut = (*theCoupleTable->GetEnergyCutsVector(1))[j];
254 
255  //create physics vectors then fill it (same parameters as lambda vector)
256  G4PhysicsVector * ptrVectorA = LambdaPhysicsVector(couple,cut);
257  G4PhysicsVector * ptrVectorB = LambdaPhysicsVector(couple,cut);
258  size_t bins = ptrVectorA->GetVectorLength();
259 
260  for (size_t i = 0 ; i < bins ; i++ ) {
261  G4double lowEdgeEnergy = ptrVectorA->Energy(i);
262  G4double tasm=0.;
263  G4double asym = ComputeAsymmetry(lowEdgeEnergy, couple, part, cut, tasm);
264  ptrVectorA->PutValue(i,asym);
265  ptrVectorB->PutValue(i,tasm);
266  }
267  theAsymmetryTable->insertAt( j , ptrVectorA ) ;
268  theTransverseAsymmetryTable->insertAt( j , ptrVectorB ) ;
269  }
270 
271 }
272 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
273 
275  const G4MaterialCutsCouple* couple,
276  const G4ParticleDefinition& aParticle,
277  G4double cut,
278  G4double & tAsymmetry)
279 {
280  G4double lAsymmetry = 0.0;
281  tAsymmetry = 0.0;
282  if (isElectron) {lAsymmetry = tAsymmetry = -1.0;}
283 
284  // calculate polarized cross section
288  G4double sigma2=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
289 
290  // calculate transversely polarized cross section
294  G4double sigma3=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
295 
296  // calculate unpolarized cross section
300  G4double sigma0=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
301  // determine assymmetries
302  if (sigma0>0.) {
303  lAsymmetry=sigma2/sigma0-1.;
304  tAsymmetry=sigma3/sigma0-1.;
305  }
306  if (std::fabs(lAsymmetry)>1.) {
307  G4cout<<" energy="<<energy<<"\n";
308  G4cout<<"WARNING lAsymmetry= "<<lAsymmetry<<" ("<<std::fabs(lAsymmetry)-1.<<")\n";
309  }
310  if (std::fabs(tAsymmetry)>1.) {
311  G4cout<<" energy="<<energy<<"\n";
312  G4cout<<"WARNING tAsymmetry= "<<tAsymmetry<<" ("<<std::fabs(tAsymmetry)-1.<<")\n";
313  }
314 // else {
315 // G4cout<<" tAsymmetry= "<<tAsymmetry<<" ("<<std::fabs(tAsymmetry)-1.<<")\n";
316 // }
317  return lAsymmetry;
318 }
319 
320 
const G4ThreeVector & GetPolarization() const
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
G4double MaxKinEnergy() const
G4bool isElectron(G4int ityp)
G4int verboseLevel
Definition: G4VProcess.hh:368
G4ePolarizedIonisation(const G4String &name="pol-eIoni")
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
void BuildPhysicsTable(const G4ParticleDefinition &)
const G4DynamicParticle * GetDynamicParticle() const
G4double ComputeAsymmetry(G4double energy, const G4MaterialCutsCouple *couple, const G4ParticleDefinition &particle, G4double cut, G4double &tasm)
G4String name
Definition: TRTMaterials.hh:40
static G4PolarizationManager * GetInstance()
G4double GetMeanFreePath(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
size_t GetVectorLength() const
void SetFluctModel(G4VEmFluctuationModel *)
G4PhysicsVector * LambdaPhysicsVector(const G4MaterialCutsCouple *, G4double cut)
G4PolarizedMollerBhabhaModel * emModel
G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
G4bool IsZero() const
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:707
G4VEmFluctuationModel * FluctModel()
G4GLOB_DLL std::ostream G4cout
virtual void InitialiseEnergyLossProcess(const G4ParticleDefinition *, const G4ParticleDefinition *)
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:483
void PutValue(size_t index, G4double theValue)
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
void AddEmModel(G4int, G4VEmModel *, G4VEmFluctuationModel *fluc=0, const G4Region *region=0)
G4double MinKinEnergy() const
G4double Energy(size_t index) const
void SetSecondaryParticle(const G4ParticleDefinition *p)
virtual void BuildPhysicsTable(const G4ParticleDefinition &)
static G4ProductionCutsTable * GetProductionCutsTable()
static G4Positron * Positron()
Definition: G4Positron.cc:94
G4LogicalVolume * GetLogicalVolume() const
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
G4double energy(const ThreeVector &p, const G4double m)
G4double GetMeanFreePath(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
static G4EmParameters * Instance()
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 *)
#define G4endl
Definition: G4ios.hh:61
G4ParticleDefinition * theElectron
double G4double
Definition: G4Types.hh:76
void SetBeamPolarization(const G4ThreeVector &pBeam)
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:714
G4ForceCondition
G4ThreeVector G4ParticleMomentum
#define DBL_MAX
Definition: templates.hh:83
static G4ThreeVector GetParticleFrameX(const G4ThreeVector &)
const G4ThreeVector & GetVolumePolarization(G4LogicalVolume *lVol) const
G4VEmFluctuationModel * flucModel
void clearAndDestroy()
G4PhysicsTable * theTransverseAsymmetryTable
size_t CurrentMaterialCutsCoupleIndex() const