Geant4  9.6.p02
 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 69847 2013-05-16 09:36:18Z 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 
53 #include "G4ePolarizedIonisation.hh"
54 #include "G4Electron.hh"
55 #include "G4UniversalFluctuation.hh"
56 #include "G4BohrFluctuations.hh"
57 #include "G4UnitsTable.hh"
58 
59 #include "G4PolarizedMollerBhabhaModel.hh"
60 #include "G4ProductionCutsTable.hh"
61 #include "G4PolarizationManager.hh"
62 #include "G4PolarizationHelper.hh"
63 #include "G4StokesVector.hh"
64 
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66 
67 G4ePolarizedIonisation::G4ePolarizedIonisation(const G4String& name)
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();
83  SetProcessSubType(fIonisation);
84  flucModel = 0;
85  emModel = 0;
86 }
87 
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
89 
90 G4ePolarizedIonisation::~G4ePolarizedIonisation()
91 {
92  if (theAsymmetryTable) {
93  theAsymmetryTable->clearAndDestroy();
94  delete theAsymmetryTable;
95  }
96  if (theTransverseAsymmetryTable) {
97  theTransverseAsymmetryTable->clearAndDestroy();
98  delete theTransverseAsymmetryTable;
99  }
100 }
101 
102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
103 
104 void G4ePolarizedIonisation::InitialiseEnergyLossProcess(
105  const G4ParticleDefinition* part,
106  const G4ParticleDefinition* /*part2*/)
107 {
108  if(!isInitialised) {
109 
110  if(part == G4Positron::Positron()) isElectron = false;
111  SetSecondaryParticle(theElectron);
112 
113 
114 
115  flucModel = new G4UniversalFluctuation();
116  //flucModel = new G4BohrFluctuations();
117 
118  // G4VEmModel* em = new G4MollerBhabhaModel();
119  emModel = new G4PolarizedMollerBhabhaModel;
120  emModel->SetLowEnergyLimit(MinKinEnergy());
121  emModel->SetHighEnergyLimit(MaxKinEnergy());
122  AddEmModel(1, emModel, flucModel);
123 
124  isInitialised = true;
125  }
126 }
127 
128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
129 
130 void G4ePolarizedIonisation::PrintInfo()
131 {
132  G4cout << " Delta cross sections from Moller+Bhabha, "
133  << "good description from 1 KeV to 100 GeV."
134  << G4endl;
135 }
136 
137 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
138 
139 G4double G4ePolarizedIonisation::GetMeanFreePath(const G4Track& track,
140  G4double step,
141  G4ForceCondition* cond)
142 {
143  // *** get unploarised mean free path from lambda table ***
144  G4double mfp = G4VEnergyLossProcess::GetMeanFreePath(track, step, cond);
145 
146 
147  // *** get asymmetry, if target is polarized ***
148  G4VPhysicalVolume* aPVolume = track.GetVolume();
149  G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
150 
151  G4PolarizationManager * polarizationManger = G4PolarizationManager::GetInstance();
152  G4bool volumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
153  const G4StokesVector ePolarization = track.GetPolarization();
154 
155  if (mfp != DBL_MAX && volumeIsPolarized && !ePolarization.IsZero()) {
156  const G4DynamicParticle* aDynamicElectron = track.GetDynamicParticle();
157  G4double eEnergy = aDynamicElectron->GetKineticEnergy();
158  const G4ParticleMomentum eDirection0 = aDynamicElectron->GetMomentumDirection();
159 
160  G4StokesVector volumePolarization = polarizationManger->GetVolumePolarization(aLVolume);
161 
162  G4bool isOutRange;
163  size_t idx = CurrentMaterialCutsCoupleIndex();
164  G4double lAsymmetry = (*theAsymmetryTable)(idx)->
165  GetValue(eEnergy, isOutRange);
166  G4double tAsymmetry = (*theTransverseAsymmetryTable)(idx)->
167  GetValue(eEnergy, isOutRange);
168 
169  // calculate longitudinal spin component
170  G4double polZZ = ePolarization.z()*
171  volumePolarization*eDirection0;
172  // calculate transvers spin components
173  G4double polXX = ePolarization.x()*
174  volumePolarization*G4PolarizationHelper::GetParticleFrameX(eDirection0);
175  G4double polYY = ePolarization.y()*
176  volumePolarization*G4PolarizationHelper::GetParticleFrameY(eDirection0);
177 
178 
179  G4double impact = 1. + polZZ*lAsymmetry + (polXX + polYY)*tAsymmetry;
180  // determine polarization dependent mean free path
181  mfp /= impact;
182  if (mfp <=0.) {
183  G4cout <<"PV impact ( "<<polXX<<" , "<<polYY<<" , "<<polZZ<<" )"<<G4endl;
184  G4cout << " impact on MFP is "<< impact <<G4endl;
185  G4cout<<" lAsymmetry= "<<lAsymmetry<<" ("<<std::fabs(lAsymmetry)-1.<<")\n";
186  G4cout<<" tAsymmetry= "<<tAsymmetry<<" ("<<std::fabs(tAsymmetry)-1.<<")\n";
187  }
188  }
189 
190  return mfp;
191 }
192 
193 G4double G4ePolarizedIonisation::PostStepGetPhysicalInteractionLength(const G4Track& track,
194  G4double step,
195  G4ForceCondition* cond)
196 {
197  // *** get unploarised mean free path from lambda table ***
198  G4double mfp = G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength(track, step, cond);
199 
200 
201  // *** get asymmetry, if target is polarized ***
202  G4VPhysicalVolume* aPVolume = track.GetVolume();
203  G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
204 
205  G4PolarizationManager * polarizationManger = G4PolarizationManager::GetInstance();
206  G4bool volumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
207  const G4StokesVector ePolarization = track.GetPolarization();
208 
209  if (mfp != DBL_MAX && volumeIsPolarized && !ePolarization.IsZero()) {
210  const G4DynamicParticle* aDynamicElectron = track.GetDynamicParticle();
211  G4double eEnergy = aDynamicElectron->GetKineticEnergy();
212  const G4ParticleMomentum eDirection0 = aDynamicElectron->GetMomentumDirection();
213 
214  G4StokesVector volumePolarization = polarizationManger->GetVolumePolarization(aLVolume);
215 
216  size_t idx = CurrentMaterialCutsCoupleIndex();
217  G4double lAsymmetry = (*theAsymmetryTable)(idx)->Value(eEnergy);
218  G4double tAsymmetry = (*theTransverseAsymmetryTable)(idx)->Value(eEnergy);
219 
220  // calculate longitudinal spin component
221  G4double polZZ = ePolarization.z()*
222  volumePolarization*eDirection0;
223  // calculate transvers spin components
224  G4double polXX = ePolarization.x()*
225  volumePolarization*G4PolarizationHelper::GetParticleFrameX(eDirection0);
226  G4double polYY = ePolarization.y()*
227  volumePolarization*G4PolarizationHelper::GetParticleFrameY(eDirection0);
228 
229 
230  G4double impact = 1. + polZZ*lAsymmetry + (polXX + polYY)*tAsymmetry;
231  // determine polarization dependent mean free path
232  mfp /= impact;
233  if (mfp <=0.) {
234  G4cout <<"PV impact ( "<<polXX<<" , "<<polYY<<" , "<<polZZ<<" )"<<G4endl;
235  G4cout << " impact on MFP is "<< impact <<G4endl;
236  G4cout<<" lAsymmetry= "<<lAsymmetry<<" ("<<std::fabs(lAsymmetry)-1.<<")\n";
237  G4cout<<" tAsymmetry= "<<tAsymmetry<<" ("<<std::fabs(tAsymmetry)-1.<<")\n";
238  }
239  }
240 
241  return mfp;
242 }
243 
244 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
245 void G4ePolarizedIonisation::BuildPhysicsTable(const G4ParticleDefinition& part)
246 {
247  // *** build DEDX and (unpolarized) cross section tables
248  G4VEnergyLossProcess::BuildPhysicsTable(part);
249  // G4PhysicsTable* pt =
250  // BuildDEDXTable();
251 
252 
253  // *** build asymmetry-table
254  if (theAsymmetryTable) {
255  theAsymmetryTable->clearAndDestroy(); delete theAsymmetryTable;}
256  if (theTransverseAsymmetryTable) {
257  theTransverseAsymmetryTable->clearAndDestroy(); delete theTransverseAsymmetryTable;}
258 
259  const G4ProductionCutsTable* theCoupleTable=
260  G4ProductionCutsTable::GetProductionCutsTable();
261  size_t numOfCouples = theCoupleTable->GetTableSize();
262 
263  theAsymmetryTable = new G4PhysicsTable(numOfCouples);
264  theTransverseAsymmetryTable = new G4PhysicsTable(numOfCouples);
265 
266  for (size_t j=0 ; j < numOfCouples; j++ ) {
267  // get cut value
268  const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
269 
270  G4double cut = (*theCoupleTable->GetEnergyCutsVector(1))[j];
271 
272  //create physics vectors then fill it (same parameters as lambda vector)
273  G4PhysicsVector * ptrVectorA = LambdaPhysicsVector(couple,cut);
274  G4PhysicsVector * ptrVectorB = LambdaPhysicsVector(couple,cut);
275  size_t bins = ptrVectorA->GetVectorLength();
276 
277  for (size_t i = 0 ; i < bins ; i++ ) {
278  G4double lowEdgeEnergy = ptrVectorA->Energy(i);
279  G4double tasm=0.;
280  G4double asym = ComputeAsymmetry(lowEdgeEnergy, couple, part, cut, tasm);
281  ptrVectorA->PutValue(i,asym);
282  ptrVectorB->PutValue(i,tasm);
283  }
284  theAsymmetryTable->insertAt( j , ptrVectorA ) ;
285  theTransverseAsymmetryTable->insertAt( j , ptrVectorB ) ;
286  }
287 
288 }
289 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
290 
291 G4double G4ePolarizedIonisation::ComputeAsymmetry(G4double energy,
292  const G4MaterialCutsCouple* couple,
293  const G4ParticleDefinition& aParticle,
294  G4double cut,
295  G4double & tAsymmetry)
296 {
297  G4double lAsymmetry = 0.0;
298  tAsymmetry = 0.0;
299  if (isElectron) {lAsymmetry = tAsymmetry = -1.0;}
300 
301  // calculate polarized cross section
302  theTargetPolarization=G4ThreeVector(0.,0.,1.);
303  emModel->SetTargetPolarization(theTargetPolarization);
304  emModel->SetBeamPolarization(theTargetPolarization);
305  G4double sigma2=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
306 
307  // calculate transversely polarized cross section
308  theTargetPolarization=G4ThreeVector(1.,0.,0.);
309  emModel->SetTargetPolarization(theTargetPolarization);
310  emModel->SetBeamPolarization(theTargetPolarization);
311  G4double sigma3=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
312 
313  // calculate unpolarized cross section
314  theTargetPolarization=G4ThreeVector();
315  emModel->SetTargetPolarization(theTargetPolarization);
316  emModel->SetBeamPolarization(theTargetPolarization);
317  G4double sigma0=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
318  // determine assymmetries
319  if (sigma0>0.) {
320  lAsymmetry=sigma2/sigma0-1.;
321  tAsymmetry=sigma3/sigma0-1.;
322  }
323  if (std::fabs(lAsymmetry)>1.) {
324  G4cout<<" energy="<<energy<<"\n";
325  G4cout<<"WARNING lAsymmetry= "<<lAsymmetry<<" ("<<std::fabs(lAsymmetry)-1.<<")\n";
326  }
327  if (std::fabs(tAsymmetry)>1.) {
328  G4cout<<" energy="<<energy<<"\n";
329  G4cout<<"WARNING tAsymmetry= "<<tAsymmetry<<" ("<<std::fabs(tAsymmetry)-1.<<")\n";
330  }
331 // else {
332 // G4cout<<" tAsymmetry= "<<tAsymmetry<<" ("<<std::fabs(tAsymmetry)-1.<<")\n";
333 // }
334  return lAsymmetry;
335 }
336 
337 
Generated on Sat May 25 2013 14:33:36 for Geant4 by   doxygen 1.8.4