Geant4_10
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 76472 2013-11-11 10:34:07Z 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 
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66 
68  : G4VEnergyLossProcess(name),
69  theElectron(G4Electron::Electron()),
70  isElectron(true),
71  isInitialised(false),
72  theAsymmetryTable(NULL),
73  theTransverseAsymmetryTable(NULL)
74 {
75  verboseLevel=0;
76  // SetDEDXBinning(120);
77  // SetLambdaBinning(120);
78  // numBinAsymmetryTable=78;
79 
80  // SetMinKinEnergy(0.1*keV);
81  // SetMaxKinEnergy(100.0*TeV);
82  // PrintInfoDefinition();
84  flucModel = 0;
85  emModel = 0;
86 }
87 
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
89 
91 {
92  delete theAsymmetryTable;
93  delete theTransverseAsymmetryTable;
94 }
95 
96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
97 
100  const G4ParticleDefinition* /*part2*/)
101 {
102  if(!isInitialised) {
103 
104  if(part == G4Positron::Positron()) isElectron = false;
105  SetSecondaryParticle(theElectron);
106 
107 
108 
109  flucModel = new G4UniversalFluctuation();
110  //flucModel = new G4BohrFluctuations();
111 
112  // G4VEmModel* em = new G4MollerBhabhaModel();
113  emModel = new G4PolarizedMollerBhabhaModel;
114  emModel->SetLowEnergyLimit(MinKinEnergy());
115  emModel->SetHighEnergyLimit(MaxKinEnergy());
116  AddEmModel(1, emModel, flucModel);
117 
118  isInitialised = true;
119  }
120 }
121 
122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
123 
125 {
126  G4cout << " Delta cross sections from Moller+Bhabha, "
127  << "good description from 1 KeV to 100 GeV."
128  << G4endl;
129 }
130 
131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
132 
134  G4double step,
135  G4ForceCondition* cond)
136 {
137  // *** get unploarised mean free path from lambda table ***
138  G4double mfp = G4VEnergyLossProcess::GetMeanFreePath(track, step, cond);
139 
140 
141  // *** get asymmetry, if target is polarized ***
142  G4VPhysicalVolume* aPVolume = track.GetVolume();
143  G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
144 
146  G4bool volumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
147  const G4StokesVector ePolarization = track.GetPolarization();
148 
149  if (mfp != DBL_MAX && volumeIsPolarized && !ePolarization.IsZero()) {
150  const G4DynamicParticle* aDynamicElectron = track.GetDynamicParticle();
151  G4double eEnergy = aDynamicElectron->GetKineticEnergy();
152  const G4ParticleMomentum eDirection0 = aDynamicElectron->GetMomentumDirection();
153 
154  G4StokesVector volumePolarization = polarizationManger->GetVolumePolarization(aLVolume);
155 
156  G4bool isOutRange;
157  size_t idx = CurrentMaterialCutsCoupleIndex();
158  G4double lAsymmetry = (*theAsymmetryTable)(idx)->
159  GetValue(eEnergy, isOutRange);
160  G4double tAsymmetry = (*theTransverseAsymmetryTable)(idx)->
161  GetValue(eEnergy, isOutRange);
162 
163  // calculate longitudinal spin component
164  G4double polZZ = ePolarization.z()*
165  volumePolarization*eDirection0;
166  // calculate transvers spin components
167  G4double polXX = ePolarization.x()*
168  volumePolarization*G4PolarizationHelper::GetParticleFrameX(eDirection0);
169  G4double polYY = ePolarization.y()*
170  volumePolarization*G4PolarizationHelper::GetParticleFrameY(eDirection0);
171 
172 
173  G4double impact = 1. + polZZ*lAsymmetry + (polXX + polYY)*tAsymmetry;
174  // determine polarization dependent mean free path
175  mfp /= impact;
176  if (mfp <=0.) {
177  G4cout <<"PV impact ( "<<polXX<<" , "<<polYY<<" , "<<polZZ<<" )"<<G4endl;
178  G4cout << " impact on MFP is "<< impact <<G4endl;
179  G4cout<<" lAsymmetry= "<<lAsymmetry<<" ("<<std::fabs(lAsymmetry)-1.<<")\n";
180  G4cout<<" tAsymmetry= "<<tAsymmetry<<" ("<<std::fabs(tAsymmetry)-1.<<")\n";
181  }
182  }
183 
184  return mfp;
185 }
186 
188  G4double step,
189  G4ForceCondition* cond)
190 {
191  // *** get unploarised mean free path from lambda table ***
193 
194 
195  // *** get asymmetry, if target is polarized ***
196  G4VPhysicalVolume* aPVolume = track.GetVolume();
197  G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
198 
200  G4bool volumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
201  const G4StokesVector ePolarization = track.GetPolarization();
202 
203  if (mfp != DBL_MAX && volumeIsPolarized && !ePolarization.IsZero()) {
204  const G4DynamicParticle* aDynamicElectron = track.GetDynamicParticle();
205  G4double eEnergy = aDynamicElectron->GetKineticEnergy();
206  const G4ParticleMomentum eDirection0 = aDynamicElectron->GetMomentumDirection();
207 
208  G4StokesVector volumePolarization = polarizationManger->GetVolumePolarization(aLVolume);
209 
210  size_t idx = CurrentMaterialCutsCoupleIndex();
211  G4double lAsymmetry = (*theAsymmetryTable)(idx)->Value(eEnergy);
212  G4double tAsymmetry = (*theTransverseAsymmetryTable)(idx)->Value(eEnergy);
213 
214  // calculate longitudinal spin component
215  G4double polZZ = ePolarization.z()*
216  volumePolarization*eDirection0;
217  // calculate transvers spin components
218  G4double polXX = ePolarization.x()*
219  volumePolarization*G4PolarizationHelper::GetParticleFrameX(eDirection0);
220  G4double polYY = ePolarization.y()*
221  volumePolarization*G4PolarizationHelper::GetParticleFrameY(eDirection0);
222 
223 
224  G4double impact = 1. + polZZ*lAsymmetry + (polXX + polYY)*tAsymmetry;
225  // determine polarization dependent mean free path
226  mfp /= impact;
227  if (mfp <=0.) {
228  G4cout <<"PV impact ( "<<polXX<<" , "<<polYY<<" , "<<polZZ<<" )"<<G4endl;
229  G4cout << " impact on MFP is "<< impact <<G4endl;
230  G4cout<<" lAsymmetry= "<<lAsymmetry<<" ("<<std::fabs(lAsymmetry)-1.<<")\n";
231  G4cout<<" tAsymmetry= "<<tAsymmetry<<" ("<<std::fabs(tAsymmetry)-1.<<")\n";
232  }
233  }
234 
235  return mfp;
236 }
237 
238 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
240 {
241  // *** build DEDX and (unpolarized) cross section tables
243  // G4PhysicsTable* pt =
244  // BuildDEDXTable();
245 
246 
247  // *** build asymmetry-table
248  if (theAsymmetryTable) {
249  theAsymmetryTable->clearAndDestroy(); delete theAsymmetryTable;}
250  if (theTransverseAsymmetryTable) {
251  theTransverseAsymmetryTable->clearAndDestroy(); delete theTransverseAsymmetryTable;}
252 
253  const G4ProductionCutsTable* theCoupleTable=
255  size_t numOfCouples = theCoupleTable->GetTableSize();
256 
257  theAsymmetryTable = new G4PhysicsTable(numOfCouples);
258  theTransverseAsymmetryTable = new G4PhysicsTable(numOfCouples);
259 
260  for (size_t j=0 ; j < numOfCouples; j++ ) {
261  // get cut value
262  const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
263 
264  G4double cut = (*theCoupleTable->GetEnergyCutsVector(1))[j];
265 
266  //create physics vectors then fill it (same parameters as lambda vector)
267  G4PhysicsVector * ptrVectorA = LambdaPhysicsVector(couple,cut);
268  G4PhysicsVector * ptrVectorB = LambdaPhysicsVector(couple,cut);
269  size_t bins = ptrVectorA->GetVectorLength();
270 
271  for (size_t i = 0 ; i < bins ; i++ ) {
272  G4double lowEdgeEnergy = ptrVectorA->Energy(i);
273  G4double tasm=0.;
274  G4double asym = ComputeAsymmetry(lowEdgeEnergy, couple, part, cut, tasm);
275  ptrVectorA->PutValue(i,asym);
276  ptrVectorB->PutValue(i,tasm);
277  }
278  theAsymmetryTable->insertAt( j , ptrVectorA ) ;
279  theTransverseAsymmetryTable->insertAt( j , ptrVectorB ) ;
280  }
281 
282 }
283 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
284 
286  const G4MaterialCutsCouple* couple,
287  const G4ParticleDefinition& aParticle,
288  G4double cut,
289  G4double & tAsymmetry)
290 {
291  G4double lAsymmetry = 0.0;
292  tAsymmetry = 0.0;
293  if (isElectron) {lAsymmetry = tAsymmetry = -1.0;}
294 
295  // calculate polarized cross section
296  theTargetPolarization=G4ThreeVector(0.,0.,1.);
297  emModel->SetTargetPolarization(theTargetPolarization);
298  emModel->SetBeamPolarization(theTargetPolarization);
299  G4double sigma2=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
300 
301  // calculate transversely polarized cross section
302  theTargetPolarization=G4ThreeVector(1.,0.,0.);
303  emModel->SetTargetPolarization(theTargetPolarization);
304  emModel->SetBeamPolarization(theTargetPolarization);
305  G4double sigma3=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
306 
307  // calculate unpolarized cross section
308  theTargetPolarization=G4ThreeVector();
309  emModel->SetTargetPolarization(theTargetPolarization);
310  emModel->SetBeamPolarization(theTargetPolarization);
311  G4double sigma0=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
312  // determine assymmetries
313  if (sigma0>0.) {
314  lAsymmetry=sigma2/sigma0-1.;
315  tAsymmetry=sigma3/sigma0-1.;
316  }
317  if (std::fabs(lAsymmetry)>1.) {
318  G4cout<<" energy="<<energy<<"\n";
319  G4cout<<"WARNING lAsymmetry= "<<lAsymmetry<<" ("<<std::fabs(lAsymmetry)-1.<<")\n";
320  }
321  if (std::fabs(tAsymmetry)>1.) {
322  G4cout<<" energy="<<energy<<"\n";
323  G4cout<<"WARNING tAsymmetry= "<<tAsymmetry<<" ("<<std::fabs(tAsymmetry)-1.<<")\n";
324  }
325 // else {
326 // G4cout<<" tAsymmetry= "<<tAsymmetry<<" ("<<std::fabs(tAsymmetry)-1.<<")\n";
327 // }
328  return lAsymmetry;
329 }
330 
331 
const G4ThreeVector & GetPolarization() const
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
G4bool isElectron(G4int ityp)
G4int verboseLevel
Definition: G4VProcess.hh:368
G4ePolarizedIonisation(const G4String &name="pol-eIoni")
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
double x() const
void BuildPhysicsTable(const G4ParticleDefinition &)
const G4DynamicParticle * GetDynamicParticle() const
G4double ComputeAsymmetry(G4double energy, const G4MaterialCutsCouple *couple, const G4ParticleDefinition &particle, G4double cut, G4double &tasm)
static G4PolarizationManager * GetInstance()
G4double GetMeanFreePath(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
const XML_Char * name
Definition: expat.h:151
G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
size_t GetVectorLength() const
G4PhysicsVector * LambdaPhysicsVector(const G4MaterialCutsCouple *, G4double cut)
G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
double z() const
G4bool IsZero() const
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:683
G4GLOB_DLL std::ostream G4cout
double energy
Definition: plottest35.C:25
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:467
TString part[npart]
Definition: Style.C:32
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 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 GetMeanFreePath(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
double y() const
bool IsPolarized(G4LogicalVolume *lVol) const
void SetTargetPolarization(const G4ThreeVector &pTarget)
G4double MinKinEnergy() const
G4VPhysicalVolume * GetVolume() const
G4double MaxKinEnergy() const
static G4ThreeVector GetParticleFrameY(const G4ThreeVector &)
void insertAt(size_t, G4PhysicsVector *)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
void SetBeamPolarization(const G4ThreeVector &pBeam)
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:690
G4ForceCondition
#define DBL_MAX
Definition: templates.hh:83
static G4ThreeVector GetParticleFrameX(const G4ThreeVector &)
const G4ThreeVector & GetVolumePolarization(G4LogicalVolume *lVol) const
void clearAndDestroy()
size_t CurrentMaterialCutsCoupleIndex() const