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