Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/polarisation/src/G4PolarizedIonisation.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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