ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4AdjointPhotoElectricModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4AdjointPhotoElectricModel.cc
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 #include "G4AdjointPhotoElectricModel.hh"
28 #include "G4AdjointCSManager.hh"
29 
30 #include "G4PhysicalConstants.hh"
31 #include "G4Integrator.hh"
32 #include "G4TrackStatus.hh"
33 #include "G4ParticleChange.hh"
34 #include "G4AdjointElectron.hh"
35 #include "G4Gamma.hh"
36 #include "G4AdjointGamma.hh"
37 
38 
40 //
41 G4AdjointPhotoElectricModel::G4AdjointPhotoElectricModel():
42  G4VEmAdjointModel("AdjointPEEffect")
43 
44 { SetUseMatrix(false);
45  SetApplyCutInRange(false);
46 
47  //Initialization
48  current_eEnergy =0.;
49  totAdjointCS=0.;
50  factorCSBiasing =1.;
51  post_step_AdjointCS =0.;
52  pre_step_AdjointCS =0.;
53  totBiasedAdjointCS =0.;
54 
55  index_element=0;
56 
57  theAdjEquivOfDirectPrimPartDef =G4AdjointGamma::AdjointGamma();
58  theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
59  theDirectPrimaryPartDef=G4Gamma::Gamma();
60  second_part_of_same_type=false;
61  theDirectPEEffectModel = new G4PEEffectFluoModel();
62 }
64 //
65 G4AdjointPhotoElectricModel::~G4AdjointPhotoElectricModel()
66 {;}
67 
69 //
70 void G4AdjointPhotoElectricModel::SampleSecondaries(const G4Track& aTrack,
71  G4bool IsScatProjToProjCase,
72  G4ParticleChange* fParticleChange)
73 { if (IsScatProjToProjCase) return ;
74 
75  //Compute the totAdjointCS vectors if not already done for the current couple and electron energy
76  //-----------------------------------------------------------------------------------------------
77  const G4MaterialCutsCouple* aCouple = aTrack.GetMaterialCutsCouple();
78  const G4DynamicParticle* aDynPart = aTrack.GetDynamicParticle() ;
79  G4double electronEnergy = aDynPart->GetKineticEnergy();
80  G4ThreeVector electronDirection= aDynPart->GetMomentumDirection() ;
81  pre_step_AdjointCS = totAdjointCS; //The last computed CS was at pre step point
82  post_step_AdjointCS = AdjointCrossSection(aCouple, electronEnergy,IsScatProjToProjCase);
83  post_step_AdjointCS = totAdjointCS;
84 
85 
86 
87 
88  //Sample element
89  //-------------
90  const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
91  size_t nelm = currentMaterial->GetNumberOfElements();
92  G4double rand_CS= G4UniformRand()*xsec[nelm-1];
93  for (index_element=0; index_element<nelm-1; index_element++){
94  if (rand_CS<xsec[index_element]) break;
95  }
96 
97  //Sample shell and binding energy
98  //-------------
99  G4int nShells = (*theElementVector)[index_element]->GetNbOfAtomicShells();
100  rand_CS= shell_prob[index_element][nShells-1]*G4UniformRand();
101  G4int i = 0;
102  for (i=0; i<nShells-1; i++){
103  if (rand_CS<shell_prob[index_element][i]) break;
104  }
105  G4double gammaEnergy= electronEnergy+(*theElementVector)[index_element]->GetAtomicShell(i);
106 
107  //Sample cos theta
108  //Copy of the G4PEEfectFluoModel cos theta sampling method ElecCosThetaDistribution.
109  //This method cannot be used directly from G4PEEfectFluoModel because it is a friend method. I should ask Vladimir to change that
110  //------------------------------------------------------------------------------------------------
111  //G4double cos_theta = theDirectPEEffectModel->ElecCosThetaDistribution(electronEnergy);
112 
113  G4double cos_theta = 1.;
114  G4double gamma = 1. + electronEnergy/electron_mass_c2;
115  if (gamma <= 5.) {
116  G4double beta = std::sqrt(gamma*gamma-1.)/gamma;
117  G4double b = 0.5*gamma*(gamma-1.)*(gamma-2);
118 
119  G4double rndm,term,greject,grejsup;
120  if (gamma < 2.) grejsup = gamma*gamma*(1.+b-beta*b);
121  else grejsup = gamma*gamma*(1.+b+beta*b);
122 
123  do { rndm = 1.-2*G4UniformRand();
124  cos_theta = (rndm+beta)/(rndm*beta+1.);
125  term = 1.-beta*cos_theta;
126  greject = (1.-cos_theta*cos_theta)*(1.+b*term)/(term*term);
127  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
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 //
168 void G4AdjointPhotoElectricModel::CorrectPostStepWeight(G4ParticleChange* fParticleChange,
169  G4double old_weight,
170  G4double adjointPrimKinEnergy,
171  G4double projectileKinEnergy ,
172  G4bool )
173 {
174  G4double new_weight=old_weight;
175 
176  G4double w_corr =G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection()/factorCSBiasing;
177  w_corr*=post_step_AdjointCS/pre_step_AdjointCS;
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 
190 G4double G4AdjointPhotoElectricModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
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();
205  for (index_element=0;index_element<nelm;index_element++){
206 
207  totAdjointCS +=AdjointCrossSectionPerAtom((*theElementVector)[index_element],electronEnergy)*theAtomNumDensityVector[index_element];
208  xsec[index_element] = totAdjointCS;
209  }
210 
211  totBiasedAdjointCS=std::min(totAdjointCS,0.01);
212 // totBiasedAdjointCS=totAdjointCS;
213  factorCSBiasing = totBiasedAdjointCS/totAdjointCS;
214  lastCS=totBiasedAdjointCS;
215 
216 
217  }
218  return totBiasedAdjointCS;
219 
220 
221 }
223 //
224 
225 G4double G4AdjointPhotoElectricModel::GetAdjointCrossSection(const G4MaterialCutsCouple* aCouple,
226  G4double electronEnergy,
227  G4bool IsScatProjToProjCase)
228 { return AdjointCrossSection(aCouple,electronEnergy,IsScatProjToProjCase);
229 }
231 //
232 
233 G4double G4AdjointPhotoElectricModel::AdjointCrossSectionPerAtom(const G4Element* anElement,G4double electronEnergy)
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;
240  G4double CS= theDirectPEEffectModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),gammaEnergy,Z,0.,0.,0.);
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 
252  CS=theDirectPEEffectModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),gammaEnergy,Z,0.,0.,0.);
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 
265 void G4AdjointPhotoElectricModel::DefineCurrentMaterialAndElectronEnergy(const G4MaterialCutsCouple* couple, G4double anEnergy)
266 { currentCouple = const_cast<G4MaterialCutsCouple*> (couple);
267  currentMaterial = const_cast<G4Material*> (couple->GetMaterial());
268  currentCoupleIndex = couple->GetIndex();
269  currentMaterialIndex = currentMaterial->GetIndex();
270  current_eEnergy = anEnergy;
271  theDirectPEEffectModel->SetCurrentCouple(couple);
272 }