Geant4  10.00.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 75591 2013-11-04 12:33:11Z 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  } while(greject < G4UniformRand()*grejsup);
129  }
130 
131  // direction of the adjoint gamma electron
132  //---------------------------------------
133 
134 
135  G4double sin_theta = std::sqrt(1.-cos_theta*cos_theta);
136  G4double Phi = twopi * G4UniformRand();
137  G4double dirx = sin_theta*std::cos(Phi),diry = sin_theta*std::sin(Phi),dirz = cos_theta;
138  G4ThreeVector adjoint_gammaDirection(dirx,diry,dirz);
139  adjoint_gammaDirection.rotateUz(electronDirection);
140 
141 
142 
143  //Weight correction
144  //-----------------------
145  CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), electronEnergy,gammaEnergy,IsScatProjToProjCase);
146 
147 
148 
149  //Create secondary and modify fParticleChange
150  //--------------------------------------------
151  G4DynamicParticle* anAdjointGamma = new G4DynamicParticle (
152  G4AdjointGamma::AdjointGamma(),adjoint_gammaDirection, gammaEnergy);
153 
154 
155 
156 
157 
158  fParticleChange->ProposeTrackStatus(fStopAndKill);
159  fParticleChange->AddSecondary(anAdjointGamma);
160 
161 
162 
163 
164 }
165 
167 //
169  G4double old_weight,
170  G4double adjointPrimKinEnergy,
171  G4double projectileKinEnergy ,
172  G4bool )
173 {
174  G4double new_weight=old_weight;
175 
178 
179 
180  new_weight*=w_corr;
181  new_weight*=projectileKinEnergy/adjointPrimKinEnergy;
182  fParticleChange->SetParentWeightByProcess(false);
183  fParticleChange->SetSecondaryWeightByProcess(false);
184  fParticleChange->ProposeParentWeight(new_weight);
185 }
186 
188 //
189 
191  G4double electronEnergy,
192  G4bool IsScatProjToProjCase)
193 {
194 
195 
196  if (IsScatProjToProjCase) return 0.;
197 
198 
199  if (aCouple !=currentCouple || current_eEnergy !=electronEnergy) {
200  totAdjointCS = 0.;
201  DefineCurrentMaterialAndElectronEnergy(aCouple, electronEnergy);
202  const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
203  const double* theAtomNumDensityVector = currentMaterial->GetVecNbOfAtomsPerVolume();
204  size_t nelm = currentMaterial->GetNumberOfElements();
206 
207  totAdjointCS +=AdjointCrossSectionPerAtom((*theElementVector)[index_element],electronEnergy)*theAtomNumDensityVector[index_element];
209  }
210 
212 // totBiasedAdjointCS=totAdjointCS;
215 
216 
217  }
218  return totBiasedAdjointCS;
219 
220 
221 }
223 //
224 
226  G4double electronEnergy,
227  G4bool IsScatProjToProjCase)
228 { return AdjointCrossSection(aCouple,electronEnergy,IsScatProjToProjCase);
229 }
231 //
232 
234 {
235  G4int nShells = anElement->GetNbOfAtomicShells();
236  G4double Z= anElement->GetZ();
237  G4int i = 0;
238  G4double B0=anElement->GetAtomicShell(0);
239  G4double gammaEnergy = electronEnergy+B0;
241  G4double adjointCS =0.;
242  if (CS >0) adjointCS += CS/gammaEnergy;
243  shell_prob[index_element][0] = adjointCS;
244  for (i=1;i<nShells;i++){
245  //G4cout<<i<<G4endl;
246  G4double Bi_= anElement->GetAtomicShell(i-1);
247  G4double Bi = anElement->GetAtomicShell(i);
248  //G4cout<<Bi_<<'\t'<<Bi<<G4endl;
249  if (electronEnergy <Bi_-Bi) {
250  gammaEnergy = electronEnergy+Bi;
251 
253  if (CS>0) adjointCS +=CS/gammaEnergy;
254  }
255  shell_prob[index_element][i] = adjointCS;
256 
257  }
258  adjointCS*=electronEnergy;
259  return adjointCS;
260 
261 }
263 //
264 
266 { currentCouple = const_cast<G4MaterialCutsCouple*> (couple);
267  currentMaterial = const_cast<G4Material*> (couple->GetMaterial());
268  currentCoupleIndex = couple->GetIndex();
270  current_eEnergy = anEnergy;
272 }
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:260
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:188
int G4int
Definition: G4Types.hh:78
void SetSecondaryWeightByProcess(G4bool)
void SetParentWeightByProcess(G4bool)
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:204
#define G4UniformRand()
Definition: Randomize.hh:87
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
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:419
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:184
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