Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/polarisation/src/G4PolarizedCompton.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 // File name:     G4PolarizedCompton
 27 //
 28 // Author:        Andreas Schaelicke
 29 //                based on code by Michel Maire / Vladimir IVANTCHENKO
 30 //
 31 // Class description
 32 //   modified version respecting media and beam polarization
 33 //   using the stokes formalism
 34 
 35 #include "G4PolarizedCompton.hh"
 36 
 37 #include "G4Electron.hh"
 38 #include "G4EmParameters.hh"
 39 #include "G4KleinNishinaCompton.hh"
 40 #include "G4PhysicsTableHelper.hh"
 41 #include "G4PolarizationManager.hh"
 42 #include "G4PolarizedComptonModel.hh"
 43 #include "G4ProductionCutsTable.hh"
 44 #include "G4StokesVector.hh"
 45 #include "G4SystemOfUnits.hh"
 46 
 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 48 G4PhysicsTable* G4PolarizedCompton::theAsymmetryTable = nullptr;
 49 
 50 G4PolarizedCompton::G4PolarizedCompton(const G4String& processName,
 51                                        G4ProcessType type)
 52   : G4VEmProcess(processName, type)
 53   , fType(10)
 54   , fBuildAsymmetryTable(true)
 55   , fUseAsymmetryTable(true)
 56   , fIsInitialised(false)
 57 {
 58   SetStartFromNullFlag(true);
 59   SetBuildTableFlag(true);
 60   SetSecondaryParticle(G4Electron::Electron());
 61   SetProcessSubType(fComptonScattering);
 62   SetMinKinEnergyPrim(1. * MeV);
 63   SetSplineFlag(true);
 64   fEmModel = nullptr;
 65 }
 66 
 67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 68 G4PolarizedCompton::~G4PolarizedCompton() { CleanTable(); }
 69 
 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 71 void G4PolarizedCompton::ProcessDescription(std::ostream& out) const
 72 {
 73   out << "Polarized model for Compton scattering.\n";
 74 
 75   G4VEmProcess::ProcessDescription(out);
 76 }
 77 
 78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 79 void G4PolarizedCompton::CleanTable()
 80 {
 81   if(theAsymmetryTable)
 82   {
 83     theAsymmetryTable->clearAndDestroy();
 84     delete theAsymmetryTable;
 85     theAsymmetryTable = nullptr;
 86   }
 87 }
 88 
 89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 90 G4bool G4PolarizedCompton::IsApplicable(const G4ParticleDefinition& p)
 91 {
 92   return (&p == G4Gamma::Gamma());
 93 }
 94 
 95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 96 void G4PolarizedCompton::InitialiseProcess(const G4ParticleDefinition*)
 97 {
 98   if(!fIsInitialised)
 99   {
100     fIsInitialised = true;
101     if(0 == fType)
102     {
103       if(nullptr == EmModel(0))
104       {
105         SetEmModel(new G4KleinNishinaCompton());
106       }
107     }
108     else
109     {
110       fEmModel = new G4PolarizedComptonModel();
111       SetEmModel(fEmModel);
112     }
113     G4EmParameters* param = G4EmParameters::Instance();
114     EmModel(0)->SetLowEnergyLimit(param->MinKinEnergy());
115     EmModel(0)->SetHighEnergyLimit(param->MaxKinEnergy());
116     AddEmModel(1, EmModel(0));
117   }
118 }
119 
120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
121 void G4PolarizedCompton::SetModel(const G4String& ss)
122 {
123   if(ss == "Klein-Nishina")
124   {
125     fType = 0;
126   }
127   if(ss == "Polarized-Compton")
128   {
129     fType = 10;
130   }
131 }
132 
133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
134 G4double G4PolarizedCompton::GetMeanFreePath(const G4Track& aTrack,
135                                              G4double previousStepSize,
136                                              G4ForceCondition* condition)
137 {
138   // *** get unploarised mean free path from lambda table ***
139   G4double mfp =
140     G4VEmProcess::GetMeanFreePath(aTrack, previousStepSize, condition);
141 
142   if(theAsymmetryTable && fUseAsymmetryTable && mfp < DBL_MAX)
143   {
144     mfp *= ComputeSaturationFactor(aTrack);
145   }
146   if(verboseLevel >= 2)
147   {
148     G4cout << "G4PolarizedCompton::MeanFreePath:  " << mfp / mm << " mm "
149            << G4endl;
150   }
151   return mfp;
152 }
153 
154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
155 G4double G4PolarizedCompton::PostStepGetPhysicalInteractionLength(
156   const G4Track& aTrack, G4double previousStepSize, G4ForceCondition* condition)
157 {
158   // save previous values
159   G4double nLength = theNumberOfInteractionLengthLeft;
160   G4double iLength = currentInteractionLength;
161 
162   // *** compute unpolarized step limit ***
163   // this changes theNumberOfInteractionLengthLeft and currentInteractionLength
164   G4double x = G4VEmProcess::PostStepGetPhysicalInteractionLength(
165     aTrack, previousStepSize, condition);
166   G4double x0      = x;
167   G4double satFact = 1.0;
168 
169   // *** add corrections on polarisation ***
170   if(theAsymmetryTable && fUseAsymmetryTable && x < DBL_MAX)
171   {
172     satFact            = ComputeSaturationFactor(aTrack);
173     G4double curLength = currentInteractionLength * satFact;
174     G4double prvLength = iLength * satFact;
175     if(nLength > 0.0)
176     {
177       theNumberOfInteractionLengthLeft =
178         std::max(nLength - previousStepSize / prvLength, 0.0);
179     }
180     x = theNumberOfInteractionLengthLeft * curLength;
181   }
182   if(verboseLevel >= 2)
183   {
184     G4cout << "G4PolarizedCompton::PostStepGPIL: " << std::setprecision(8)
185            << x / mm << " mm;" << G4endl
186            << "               unpolarized value: " << std::setprecision(8)
187            << x0 / mm << " mm." << G4endl;
188   }
189   return x;
190 }
191 
192 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
193 G4double G4PolarizedCompton::ComputeSaturationFactor(const G4Track& aTrack)
194 {
195   G4double factor = 1.0;
196 
197   // *** get asymmetry, if target is polarized ***
198   const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
199   const G4double GammaEnergy             = aDynamicGamma->GetKineticEnergy();
200   const G4StokesVector GammaPolarization =
201     G4StokesVector(aTrack.GetPolarization());
202   const G4ParticleMomentum GammaDirection0 =
203     aDynamicGamma->GetMomentumDirection();
204 
205   const G4Material* aMaterial = aTrack.GetMaterial();
206   G4VPhysicalVolume* aPVolume = aTrack.GetVolume();
207   G4LogicalVolume* aLVolume   = aPVolume->GetLogicalVolume();
208 
209   G4PolarizationManager* polarizationManager =
210     G4PolarizationManager::GetInstance();
211 
212   const G4bool VolumeIsPolarized = polarizationManager->IsPolarized(aLVolume);
213   G4StokesVector ElectronPolarization =
214     polarizationManager->GetVolumePolarization(aLVolume);
215 
216   if(VolumeIsPolarized)
217   {
218     if(verboseLevel >= 2)
219     {
220       G4cout << "G4PolarizedCompton::ComputeSaturationFactor: " << G4endl;
221       G4cout << " Mom " << GammaDirection0 << G4endl;
222       G4cout << " Polarization " << GammaPolarization << G4endl;
223       G4cout << " MaterialPol. " << ElectronPolarization << G4endl;
224       G4cout << " Phys. Volume " << aPVolume->GetName() << G4endl;
225       G4cout << " Log. Volume  " << aLVolume->GetName() << G4endl;
226       G4cout << " Material     " << aMaterial << G4endl;
227     }
228 
229     std::size_t midx               = CurrentMaterialCutsCoupleIndex();
230     const G4PhysicsVector* aVector = nullptr;
231     if(midx < theAsymmetryTable->size())
232     {
233       aVector = (*theAsymmetryTable)(midx);
234     }
235     if(aVector)
236     {
237       G4double asymmetry = aVector->Value(GammaEnergy);
238 
239       //  we have to determine angle between particle motion
240       //  and target polarisation here
241       //      circ pol * Vec(ElectronPol)*Vec(PhotonMomentum)
242       //  both vectors in global reference frame
243 
244       G4double pol        = ElectronPolarization * GammaDirection0;
245       G4double polProduct = GammaPolarization.p3() * pol;
246       factor /= (1. + polProduct * asymmetry);
247       if(verboseLevel >= 2)
248       {
249         G4cout << " Asymmetry:     " << asymmetry << G4endl;
250         G4cout << " PolProduct:    " << polProduct << G4endl;
251         G4cout << " Factor:        " << factor << G4endl;
252       }
253     }
254     else
255     {
256       G4ExceptionDescription ed;
257       ed << "Problem with asymmetry table: material index " << midx
258          << " is out of range or the table is not filled";
259       G4Exception("G4PolarizedComptonModel::ComputeSaturationFactor", "em0048",
260                   JustWarning, ed, "");
261     }
262   }
263   return factor;
264 }
265 
266 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
267 void G4PolarizedCompton::BuildPhysicsTable(const G4ParticleDefinition& part)
268 {
269   // *** build (unpolarized) cross section tables (Lambda)
270   G4VEmProcess::BuildPhysicsTable(part);
271   if(fBuildAsymmetryTable && fEmModel)
272   {
273     G4bool isMaster = true;
274     const G4PolarizedCompton* masterProcess =
275       static_cast<const G4PolarizedCompton*>(GetMasterProcess());
276     if(masterProcess && masterProcess != this)
277     {
278       isMaster = false;
279     }
280     if(isMaster)
281     {
282       BuildAsymmetryTable(part);
283     }
284   }
285 }
286 
287 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
288 void G4PolarizedCompton::BuildAsymmetryTable(const G4ParticleDefinition& part)
289 {
290   // cleanup old, initialise new table
291   CleanTable();
292   theAsymmetryTable =
293     G4PhysicsTableHelper::PreparePhysicsTable(theAsymmetryTable);
294 
295   // Access to materials
296   const G4ProductionCutsTable* theCoupleTable =
297     G4ProductionCutsTable::GetProductionCutsTable();
298   G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
299   if(!theAsymmetryTable)
300   {
301     return;
302   }
303   G4int nbins                 = LambdaBinning();
304   G4double emin               = MinKinEnergy();
305   G4double emax               = MaxKinEnergy();
306   G4PhysicsLogVector* aVector = nullptr;
307   G4PhysicsLogVector* bVector = nullptr;
308 
309   for(G4int i = 0; i < numOfCouples; ++i)
310   {
311     if(theAsymmetryTable->GetFlag(i))
312     {
313       // create physics vector and fill it
314       const G4MaterialCutsCouple* couple =
315         theCoupleTable->GetMaterialCutsCouple(i);
316       // use same parameters as for lambda
317       if(!aVector)
318       {
319         aVector = new G4PhysicsLogVector(emin, emax, nbins, true);
320         bVector = aVector;
321       }
322       else
323       {
324         bVector = new G4PhysicsLogVector(*aVector);
325       }
326 
327       for(G4int j = 0; j <= nbins; ++j)
328       {
329         G4double energy = bVector->Energy(j);
330         G4double tasm   = 0.;
331         G4double asym   = ComputeAsymmetry(energy, couple, part, 0., tasm);
332         bVector->PutValue(j, asym);
333       }
334       bVector->FillSecondDerivatives();
335       G4PhysicsTableHelper::SetPhysicsVector(theAsymmetryTable, i, bVector);
336     }
337   }
338 }
339 
340 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
341 G4double G4PolarizedCompton::ComputeAsymmetry(
342   G4double energy, const G4MaterialCutsCouple* couple,
343   const G4ParticleDefinition& aParticle, G4double cut, G4double& tAsymmetry)
344 {
345   G4double lAsymmetry = 0.0;
346   tAsymmetry          = 0;
347 
348   // calculate polarized cross section
349   G4ThreeVector thePolarization = G4ThreeVector(0., 0., 1.);
350   fEmModel->SetTargetPolarization(thePolarization);
351   fEmModel->SetBeamPolarization(thePolarization);
352   G4double sigma2 =
353     fEmModel->CrossSection(couple, &aParticle, energy, cut, energy);
354 
355   // calculate unpolarized cross section
356   thePolarization = G4ThreeVector();
357   fEmModel->SetTargetPolarization(thePolarization);
358   fEmModel->SetBeamPolarization(thePolarization);
359   G4double sigma0 =
360     fEmModel->CrossSection(couple, &aParticle, energy, cut, energy);
361 
362   // determine asymmetries
363   if(sigma0 > 0.)
364   {
365     lAsymmetry = sigma2 / sigma0 - 1.;
366   }
367   return lAsymmetry;
368 }
369