Geant4  10.02.p01
G4AdjointPhotoElectricModel.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 // $Id: G4AdjointPhotoElectricModel.cc 91870 2015-08-07 15:21:40Z gcosmo $
27 //
29 #include "G4AdjointCSManager.hh"
30 
31 #include "G4PhysicalConstants.hh"
32 #include "G4Integrator.hh"
33 #include "G4TrackStatus.hh"
34 #include "G4ParticleChange.hh"
35 #include "G4AdjointElectron.hh"
36 #include "G4Gamma.hh"
37 #include "G4AdjointGamma.hh"
38 
39 
41 //
43  G4VEmAdjointModel("AdjointPEEffect")
44 
45 { SetUseMatrix(false);
46  SetApplyCutInRange(false);
47 
48  //Initialization
49  current_eEnergy =0.;
50  totAdjointCS=0.;
51  factorCSBiasing =1.;
55 
56  index_element=0;
57 
63 }
65 //
67 {;}
68 
70 //
72  G4bool IsScatProjToProjCase,
73  G4ParticleChange* fParticleChange)
74 { if (IsScatProjToProjCase) return ;
75 
76  //Compute the totAdjointCS vectors if not already done for the current couple and electron energy
77  //-----------------------------------------------------------------------------------------------
78  const G4MaterialCutsCouple* aCouple = aTrack.GetMaterialCutsCouple();
79  const G4DynamicParticle* aDynPart = aTrack.GetDynamicParticle() ;
80  G4double electronEnergy = aDynPart->GetKineticEnergy();
81  G4ThreeVector electronDirection= aDynPart->GetMomentumDirection() ;
82  pre_step_AdjointCS = totAdjointCS; //The last computed CS was at pre step point
83  post_step_AdjointCS = AdjointCrossSection(aCouple, electronEnergy,IsScatProjToProjCase);
85 
86 
87 
88 
89  //Sample element
90  //-------------
91  const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
92  size_t nelm = currentMaterial->GetNumberOfElements();
93  G4double rand_CS= G4UniformRand()*xsec[nelm-1];
94  for (index_element=0; index_element<nelm-1; index_element++){
95  if (rand_CS<xsec[index_element]) break;
96  }
97 
98  //Sample shell and binding energy
99  //-------------
100  G4int nShells = (*theElementVector)[index_element]->GetNbOfAtomicShells();
101  rand_CS= shell_prob[index_element][nShells-1]*G4UniformRand();
102  G4int i = 0;
103  for (i=0; i<nShells-1; i++){
104  if (rand_CS<shell_prob[index_element][i]) break;
105  }
106  G4double gammaEnergy= electronEnergy+(*theElementVector)[index_element]->GetAtomicShell(i);
107 
108  //Sample cos theta
109  //Copy of the G4PEEfectFluoModel cos theta sampling method ElecCosThetaDistribution.
110  //This method cannot be used directly from G4PEEfectFluoModel because it is a friend method. I should ask Vladimir to change that
111  //------------------------------------------------------------------------------------------------
112  //G4double cos_theta = theDirectPEEffectModel->ElecCosThetaDistribution(electronEnergy);
113 
114  G4double cos_theta = 1.;
115  G4double gamma = 1. + electronEnergy/electron_mass_c2;
116  if (gamma <= 5.) {
117  G4double beta = std::sqrt(gamma*gamma-1.)/gamma;
118  G4double b = 0.5*gamma*(gamma-1.)*(gamma-2);
119 
120  G4double rndm,term,greject,grejsup;
121  if (gamma < 2.) grejsup = gamma*gamma*(1.+b-beta*b);
122  else grejsup = gamma*gamma*(1.+b+beta*b);
123 
124  do { rndm = 1.-2*G4UniformRand();
125  cos_theta = (rndm+beta)/(rndm*beta+1.);
126  term = 1.-beta*cos_theta;
127  greject = (1.-cos_theta*cos_theta)*(1.+b*term)/(term*term);
128  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
129  } while(greject < G4UniformRand()*grejsup);
130  }
131 
132  // direction of the adjoint gamma electron
133  //---------------------------------------
134 
135 
136  G4double sin_theta = std::sqrt(1.-cos_theta*cos_theta);
137  G4double Phi = twopi * G4UniformRand();
138  G4double dirx = sin_theta*std::cos(Phi),diry = sin_theta*std::sin(Phi),dirz = cos_theta;
139  G4ThreeVector adjoint_gammaDirection(dirx,diry,dirz);
140  adjoint_gammaDirection.rotateUz(electronDirection);
141 
142 
143 
144  //Weight correction
145  //-----------------------
146  CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), electronEnergy,gammaEnergy,IsScatProjToProjCase);
147 
148 
149 
150  //Create secondary and modify fParticleChange
151  //--------------------------------------------
152  G4DynamicParticle* anAdjointGamma = new G4DynamicParticle (
153  G4AdjointGamma::AdjointGamma(),adjoint_gammaDirection, gammaEnergy);
154 
155 
156 
157 
158 
159  fParticleChange->ProposeTrackStatus(fStopAndKill);
160  fParticleChange->AddSecondary(anAdjointGamma);
161 
162 
163 
164 
165 }
166 
168 //
170  G4double old_weight,
171  G4double adjointPrimKinEnergy,
172  G4double projectileKinEnergy ,
173  G4bool )
174 {
175  G4double new_weight=old_weight;
176 
179 
180 
181  new_weight*=w_corr;
182  new_weight*=projectileKinEnergy/adjointPrimKinEnergy;
183  fParticleChange->SetParentWeightByProcess(false);
184  fParticleChange->SetSecondaryWeightByProcess(false);
185  fParticleChange->ProposeParentWeight(new_weight);
186 }
187 
189 //
190 
192  G4double electronEnergy,
193  G4bool IsScatProjToProjCase)
194 {
195 
196 
197  if (IsScatProjToProjCase) return 0.;
198 
199 
200  if (aCouple !=currentCouple || current_eEnergy !=electronEnergy) {
201  totAdjointCS = 0.;
202  DefineCurrentMaterialAndElectronEnergy(aCouple, electronEnergy);
203  const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
204  const double* theAtomNumDensityVector = currentMaterial->GetVecNbOfAtomsPerVolume();
205  size_t nelm = currentMaterial->GetNumberOfElements();
207 
208  totAdjointCS +=AdjointCrossSectionPerAtom((*theElementVector)[index_element],electronEnergy)*theAtomNumDensityVector[index_element];
210  }
211 
213 // totBiasedAdjointCS=totAdjointCS;
216 
217 
218  }
219  return totBiasedAdjointCS;
220 
221 
222 }
224 //
225 
227  G4double electronEnergy,
228  G4bool IsScatProjToProjCase)
229 { return AdjointCrossSection(aCouple,electronEnergy,IsScatProjToProjCase);
230 }
232 //
233 
235 {
236  G4int nShells = anElement->GetNbOfAtomicShells();
237  G4double Z= anElement->GetZ();
238  G4int i = 0;
239  G4double B0=anElement->GetAtomicShell(0);
240  G4double gammaEnergy = electronEnergy+B0;
242  G4double adjointCS =0.;
243  if (CS >0) adjointCS += CS/gammaEnergy;
244  shell_prob[index_element][0] = adjointCS;
245  for (i=1;i<nShells;i++){
246  //G4cout<<i<<G4endl;
247  G4double Bi_= anElement->GetAtomicShell(i-1);
248  G4double Bi = anElement->GetAtomicShell(i);
249  //G4cout<<Bi_<<'\t'<<Bi<<G4endl;
250  if (electronEnergy <Bi_-Bi) {
251  gammaEnergy = electronEnergy+Bi;
252 
254  if (CS>0) adjointCS +=CS/gammaEnergy;
255  }
256  shell_prob[index_element][i] = adjointCS;
257 
258  }
259  adjointCS*=electronEnergy;
260  return adjointCS;
261 
262 }
264 //
265 
267 { currentCouple = const_cast<G4MaterialCutsCouple*> (couple);
268  currentMaterial = const_cast<G4Material*> (couple->GetMaterial());
269  currentCoupleIndex = couple->GetIndex();
271  current_eEnergy = anEnergy;
273 }
static G4AdjointGamma * AdjointGamma()
void DefineCurrentMaterialAndElectronEnergy(const G4MaterialCutsCouple *aCouple, G4double eEnergy)
G4int GetNbOfAtomicShells() const
Definition: G4Element.hh:146
virtual void CorrectPostStepWeight(G4ParticleChange *fParticleChange, G4double old_weight, G4double adjointPrimKinEnergy, G4double projectileKinEnergy, G4bool IsScatProjToProjCase)
std::vector< G4Element * > G4ElementVector
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double GetPostStepWeightCorrection()
const G4DynamicParticle * GetDynamicParticle() const
G4ParticleDefinition * theDirectPrimaryPartDef
size_t GetIndex() const
Definition: G4Material.hh:262
void ProposeParentWeight(G4double finalWeight)
G4double GetZ() const
Definition: G4Element.hh:131
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
static G4AdjointElectron * AdjointElectron()
virtual G4double GetAdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
G4Material * currentMaterial
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
void SetSecondaryWeightByProcess(G4bool)
void SetParentWeightByProcess(G4bool)
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:206
#define G4UniformRand()
Definition: Randomize.hh:97
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
static const double twopi
Definition: G4SIunits.hh:75
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A, G4double, G4double)
G4double AdjointCrossSectionPerAtom(const G4Element *anElement, G4double electronEnergy)
void SetUseMatrix(G4bool aBool)
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4MaterialCutsCouple * currentCouple
virtual void SampleSecondaries(const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
void SetCurrentCouple(const G4MaterialCutsCouple *)
Definition: G4VEmModel.hh:445
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
void AddSecondary(G4Track *aSecondary)
G4double GetWeight() const
G4ParticleDefinition * theAdjEquivOfDirectSecondPartDef
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
G4double GetAtomicShell(G4int index) const
Definition: G4Element.cc:365
void SetApplyCutInRange(G4bool aBool)
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
static G4AdjointCSManager * GetAdjointCSManager()
const G4Material * GetMaterial() const