source: trunk/source/processes/electromagnetic/polarisation/src/G4ePolarizedIonisation.cc @ 1350

Last change on this file since 1350 was 1337, checked in by garnier, 14 years ago

tag geant4.9.4 beta 1 + modifs locales
File size: 13.0 KB
Line 
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,v 1.8 2010/06/16 11:20:54 schaelic Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-01 $
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name:     G4ePolarizedIonisation
35//
36// Author:        A.Schaelicke on base of Vladimir Ivanchenko code
37//
38// Creation date: 10.11.2005
39//
40// Modifications:
41//
42// 10-11-05, include polarization description (A.Schaelicke)
43// , create asymmetry table and determine interactionlength
44// , update polarized differential cross section
45//
46// 20-08-06, modified interface (A.Schaelicke)
47// 11-06-07, add PostStepGetPhysicalInteractionLength (A.Schalicke)
48//
49// Class Description:
50//
51//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
52//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
53
54#include "G4ePolarizedIonisation.hh"
55#include "G4Electron.hh"
56#include "G4UniversalFluctuation.hh"
57#include "G4BohrFluctuations.hh"
58#include "G4UnitsTable.hh"
59
60#include "G4PolarizedMollerBhabhaModel.hh"
61#include "G4ProductionCutsTable.hh"
62#include "G4PolarizationManager.hh"
63#include "G4PolarizationHelper.hh"
64#include "G4StokesVector.hh"
65
66//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
67
68G4ePolarizedIonisation::G4ePolarizedIonisation(const G4String& name)
69  : G4VEnergyLossProcess(name),
70    theElectron(G4Electron::Electron()),
71    isElectron(true),
72    isInitialised(false),
73    theAsymmetryTable(NULL),
74    theTransverseAsymmetryTable(NULL)
75{
76  verboseLevel=0;
77  //  SetDEDXBinning(120);
78  //  SetLambdaBinning(120);
79  //  numBinAsymmetryTable=78;
80
81  //  SetMinKinEnergy(0.1*keV);
82  //  SetMaxKinEnergy(100.0*TeV);
83  //  PrintInfoDefinition();
84  SetProcessSubType(fIonisation);
85}
86
87//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
88
89G4ePolarizedIonisation::~G4ePolarizedIonisation()
90{
91  if (theAsymmetryTable) {
92    theAsymmetryTable->clearAndDestroy();
93    delete theAsymmetryTable;
94  }
95  if (theTransverseAsymmetryTable) {
96    theTransverseAsymmetryTable->clearAndDestroy();
97    delete theTransverseAsymmetryTable;
98  }
99}
100
101//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
102
103void G4ePolarizedIonisation::InitialiseEnergyLossProcess(
104                    const G4ParticleDefinition* part,
105                    const G4ParticleDefinition* /*part2*/)
106{
107  if(!isInitialised) {
108
109    if(part == G4Positron::Positron()) isElectron = false;
110    SetSecondaryParticle(theElectron);
111
112
113
114    flucModel = new G4UniversalFluctuation();
115    //flucModel = new G4BohrFluctuations();
116
117    //    G4VEmModel* em = new G4MollerBhabhaModel();
118    emModel = new G4PolarizedMollerBhabhaModel;
119    emModel->SetLowEnergyLimit(MinKinEnergy());
120    emModel->SetHighEnergyLimit(MaxKinEnergy());
121    AddEmModel(1, emModel, flucModel);
122
123    isInitialised = true;
124  }
125}
126
127//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
128
129void G4ePolarizedIonisation::PrintInfo()
130{
131  G4cout << "      Delta cross sections from Moller+Bhabha, "
132         << "good description from 1 KeV to 100 GeV."
133         << G4endl;
134}
135
136//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
137
138G4double G4ePolarizedIonisation::GetMeanFreePath(const G4Track& track,
139                                              G4double s,
140                                              G4ForceCondition* cond)
141{
142  // *** get unploarised mean free path from lambda table ***
143  G4double mfp = G4VEnergyLossProcess::GetMeanFreePath(track, s, cond);
144
145
146  // *** get asymmetry, if target is polarized ***
147  G4VPhysicalVolume*  aPVolume  = track.GetVolume();
148  G4LogicalVolume*    aLVolume  = aPVolume->GetLogicalVolume();
149
150  G4PolarizationManager * polarizationManger = G4PolarizationManager::GetInstance();
151  const G4bool volumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
152  const G4StokesVector ePolarization = track.GetPolarization();
153
154  if (mfp != DBL_MAX &&  volumeIsPolarized && !ePolarization.IsZero()) {
155    const G4DynamicParticle* aDynamicElectron = track.GetDynamicParticle();
156    const G4double eEnergy = aDynamicElectron->GetKineticEnergy();
157    const G4ParticleMomentum eDirection0 = aDynamicElectron->GetMomentumDirection();
158
159    G4StokesVector volumePolarization = polarizationManger->GetVolumePolarization(aLVolume);
160
161    G4bool isOutRange;
162    size_t idx = CurrentMaterialCutsCoupleIndex();
163    G4double lAsymmetry = (*theAsymmetryTable)(idx)->
164                                  GetValue(eEnergy, isOutRange);
165    G4double tAsymmetry = (*theTransverseAsymmetryTable)(idx)->
166                                  GetValue(eEnergy, isOutRange);
167
168    // calculate longitudinal spin component
169    G4double polZZ = ePolarization.z()*
170                        volumePolarization*eDirection0;
171    // calculate transvers spin components
172    G4double polXX = ePolarization.x()*
173                        volumePolarization*G4PolarizationHelper::GetParticleFrameX(eDirection0);
174    G4double polYY = ePolarization.y()*
175                        volumePolarization*G4PolarizationHelper::GetParticleFrameY(eDirection0);
176
177
178    G4double impact = 1. + polZZ*lAsymmetry + (polXX + polYY)*tAsymmetry;
179    // determine polarization dependent mean free path
180    mfp /= impact;
181    if (mfp <=0.) {
182     G4cout <<"PV impact ( "<<polXX<<" , "<<polYY<<" , "<<polZZ<<" )"<<G4endl;
183     G4cout << " impact on MFP is "<< impact <<G4endl;
184     G4cout<<" lAsymmetry= "<<lAsymmetry<<" ("<<std::fabs(lAsymmetry)-1.<<")\n";
185     G4cout<<" tAsymmetry= "<<tAsymmetry<<" ("<<std::fabs(tAsymmetry)-1.<<")\n";
186    }
187  }
188
189  return mfp;
190}
191
192G4double G4ePolarizedIonisation::PostStepGetPhysicalInteractionLength(const G4Track& track,
193                                              G4double s,
194                                              G4ForceCondition* cond)
195{
196  // *** get unploarised mean free path from lambda table ***
197  G4double mfp = G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength(track, s, cond);
198
199
200  // *** get asymmetry, if target is polarized ***
201  G4VPhysicalVolume*  aPVolume  = track.GetVolume();
202  G4LogicalVolume*    aLVolume  = aPVolume->GetLogicalVolume();
203
204  G4PolarizationManager * polarizationManger = G4PolarizationManager::GetInstance();
205  const G4bool volumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
206  const G4StokesVector ePolarization = track.GetPolarization();
207
208  if (mfp != DBL_MAX &&  volumeIsPolarized && !ePolarization.IsZero()) {
209    const G4DynamicParticle* aDynamicElectron = track.GetDynamicParticle();
210    const G4double eEnergy = aDynamicElectron->GetKineticEnergy();
211    const G4ParticleMomentum eDirection0 = aDynamicElectron->GetMomentumDirection();
212
213    G4StokesVector volumePolarization = polarizationManger->GetVolumePolarization(aLVolume);
214
215    G4bool isOutRange;
216    size_t idx = CurrentMaterialCutsCoupleIndex();
217    G4double lAsymmetry = (*theAsymmetryTable)(idx)->
218                                  GetValue(eEnergy, isOutRange);
219    G4double tAsymmetry = (*theTransverseAsymmetryTable)(idx)->
220                                  GetValue(eEnergy, isOutRange);
221
222    // calculate longitudinal spin component
223    G4double polZZ = ePolarization.z()*
224                        volumePolarization*eDirection0;
225    // calculate transvers spin components
226    G4double polXX = ePolarization.x()*
227                        volumePolarization*G4PolarizationHelper::GetParticleFrameX(eDirection0);
228    G4double polYY = ePolarization.y()*
229                        volumePolarization*G4PolarizationHelper::GetParticleFrameY(eDirection0);
230
231
232    G4double impact = 1. + polZZ*lAsymmetry + (polXX + polYY)*tAsymmetry;
233    // determine polarization dependent mean free path
234    mfp /= impact;
235    if (mfp <=0.) {
236     G4cout <<"PV impact ( "<<polXX<<" , "<<polYY<<" , "<<polZZ<<" )"<<G4endl;
237     G4cout << " impact on MFP is "<< impact <<G4endl;
238     G4cout<<" lAsymmetry= "<<lAsymmetry<<" ("<<std::fabs(lAsymmetry)-1.<<")\n";
239     G4cout<<" tAsymmetry= "<<tAsymmetry<<" ("<<std::fabs(tAsymmetry)-1.<<")\n";
240    }
241  }
242
243  return mfp;
244}
245
246//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
247void G4ePolarizedIonisation::BuildPhysicsTable(const G4ParticleDefinition& part)
248{
249  // *** build DEDX and (unpolarized) cross section tables
250  G4VEnergyLossProcess::BuildPhysicsTable(part);
251  //  G4PhysicsTable* pt =
252  //  BuildDEDXTable();
253
254
255  // *** build asymmetry-table
256  if (theAsymmetryTable) {
257    theAsymmetryTable->clearAndDestroy(); delete theAsymmetryTable;}
258  if (theTransverseAsymmetryTable) {
259    theTransverseAsymmetryTable->clearAndDestroy(); delete theTransverseAsymmetryTable;}
260
261  const G4ProductionCutsTable* theCoupleTable=
262        G4ProductionCutsTable::GetProductionCutsTable();
263  size_t numOfCouples = theCoupleTable->GetTableSize();
264
265  theAsymmetryTable = new G4PhysicsTable(numOfCouples);
266  theTransverseAsymmetryTable = new G4PhysicsTable(numOfCouples);
267
268  for (size_t j=0 ; j < numOfCouples; j++ ) {
269    // get cut value
270    const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
271
272    G4double tcutmin = emModel->MinEnergyCut(&part, couple);
273    G4double cut = (*theCoupleTable->GetEnergyCutsVector(1))[j];
274    cut = std::max(cut, tcutmin);
275
276    //create physics vectors then fill it (same parameters as lambda vector)
277    G4PhysicsVector * ptrVectorA = LambdaPhysicsVector(couple,cut);
278    G4PhysicsVector * ptrVectorB = LambdaPhysicsVector(couple,cut);
279    size_t nBins = ptrVectorA->GetVectorLength();
280
281    for (size_t i = 0 ; i < nBins ; i++ ) {
282      G4double lowEdgeEnergy = ptrVectorA->GetLowEdgeEnergy(i);
283      G4double tasm=0.;
284      G4double asym = ComputeAsymmetry(lowEdgeEnergy, couple, part, cut, tasm);
285      ptrVectorA->PutValue(i,asym);
286      ptrVectorB->PutValue(i,tasm);
287    }
288    theAsymmetryTable->insertAt( j , ptrVectorA ) ;
289    theTransverseAsymmetryTable->insertAt( j , ptrVectorB ) ;
290  }
291
292}
293//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
294
295G4double G4ePolarizedIonisation::ComputeAsymmetry(G4double energy,
296                                         const G4MaterialCutsCouple* couple,
297                                               const G4ParticleDefinition& aParticle,
298                                               G4double cut,
299                                               G4double & tAsymmetry)
300{
301  G4double lAsymmetry = 0.0;
302           tAsymmetry = 0.0;
303  if (isElectron) {lAsymmetry = tAsymmetry = -1.0;}
304
305  // calculate polarized cross section
306  theTargetPolarization=G4ThreeVector(0.,0.,1.);
307  emModel->SetTargetPolarization(theTargetPolarization);
308  emModel->SetBeamPolarization(theTargetPolarization);
309  G4double sigma2=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
310
311  // calculate transversely polarized cross section
312  theTargetPolarization=G4ThreeVector(1.,0.,0.);
313  emModel->SetTargetPolarization(theTargetPolarization);
314  emModel->SetBeamPolarization(theTargetPolarization);
315  G4double sigma3=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
316
317  // calculate unpolarized cross section
318  theTargetPolarization=G4ThreeVector();
319  emModel->SetTargetPolarization(theTargetPolarization);
320  emModel->SetBeamPolarization(theTargetPolarization);
321  G4double sigma0=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
322  // determine assymmetries
323  if (sigma0>0.) {
324    lAsymmetry=sigma2/sigma0-1.;
325    tAsymmetry=sigma3/sigma0-1.;
326  }
327  if (std::fabs(lAsymmetry)>1.) {
328    G4cout<<" energy="<<energy<<"\n";
329    G4cout<<"WARNING lAsymmetry= "<<lAsymmetry<<" ("<<std::fabs(lAsymmetry)-1.<<")\n";
330  }
331  if (std::fabs(tAsymmetry)>1.) {
332    G4cout<<" energy="<<energy<<"\n";
333    G4cout<<"WARNING tAsymmetry= "<<tAsymmetry<<" ("<<std::fabs(tAsymmetry)-1.<<")\n";
334  }
335//   else {
336//     G4cout<<"        tAsymmetry= "<<tAsymmetry<<" ("<<std::fabs(tAsymmetry)-1.<<")\n";
337//   }
338  return lAsymmetry;
339}
340
341
Note: See TracBrowser for help on using the repository browser.