Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4VEmAdjointModel.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: G4VEmAdjointModel.cc 69844 2013-05-16 09:19:33Z gcosmo $
27 //
28 #include "G4VEmAdjointModel.hh"
29 #include "G4AdjointCSManager.hh"
30 #include "G4Integrator.hh"
31 #include "G4TrackStatus.hh"
32 #include "G4ParticleChange.hh"
33 #include "G4AdjointElectron.hh"
34 #include "G4AdjointGamma.hh"
35 #include "G4AdjointPositron.hh"
36 #include "G4AdjointInterpolator.hh"
37 #include "G4PhysicsTable.hh"
38 
40 //
41 G4VEmAdjointModel::G4VEmAdjointModel(const G4String& nam):
42 name(nam)
43 // lowLimit(0.1*keV), highLimit(100.0*TeV), fluc(0), name(nam), pParticleChange(0)
44 {
45  model_index = G4AdjointCSManager::GetAdjointCSManager()->RegisterEmAdjointModel(this);
46  second_part_of_same_type =false;
47  theDirectEMModel=0;
48  mass_ratio_product=1.;
49  mass_ratio_projectile=1.;
50  currentCouple=0;
51 }
53 //
54 G4VEmAdjointModel::~G4VEmAdjointModel()
55 {;}
57 //
58 G4double G4VEmAdjointModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
59  G4double primEnergy,
60  G4bool IsScatProjToProjCase)
61 {
62  DefineCurrentMaterial(aCouple);
63  preStepEnergy=primEnergy;
64 
65  std::vector<G4double>* CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
66  if (IsScatProjToProjCase) CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
67  lastCS = G4AdjointCSManager::GetAdjointCSManager()->ComputeAdjointCS(currentMaterial,
68  this,
69  primEnergy,
70  currentTcutForDirectSecond,
71  IsScatProjToProjCase,
72  *CS_Vs_Element);
73  if (IsScatProjToProjCase) lastAdjointCSForScatProjToProjCase = lastCS;
74  else lastAdjointCSForProdToProjCase =lastCS;
75 
76 
77 
78  return lastCS;
79 
80 }
82 //
83 G4double G4VEmAdjointModel::GetAdjointCrossSection(const G4MaterialCutsCouple* aCouple,
84  G4double primEnergy,
85  G4bool IsScatProjToProjCase)
86 {
87  return AdjointCrossSection(aCouple, primEnergy,
88  IsScatProjToProjCase);
89 
90  /*
91  //To continue
92  DefineCurrentMaterial(aCouple);
93  preStepEnergy=primEnergy;
94  if (IsScatProjToProjCase){
95  G4double ekin=primEnergy*mass_ratio_projectile;
96  lastCS = G4AdjointCSManager::GetAdjointCSManager()->GetAdjointSigma(ekin, model_index,true, aCouple);
97  lastAdjointCSForScatProjToProjCase = lastCS;
98  //G4cout<<ekin<<std::endl;
99  }
100  else {
101  G4double ekin=primEnergy*mass_ratio_product;
102  lastCS = G4AdjointCSManager::GetAdjointCSManager()->GetAdjointSigma(ekin, model_index,false, aCouple);
103  lastAdjointCSForProdToProjCase = lastCS;
104  //G4cout<<ekin<<std::endl;
105  }
106  return lastCS;
107  */
108 }
110 //
111 //General implementation correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
112 G4double G4VEmAdjointModel::DiffCrossSectionPerAtomPrimToSecond(
113  G4double kinEnergyProj,
114  G4double kinEnergyProd,
115  G4double Z,
116  G4double A)
117 {
118  G4double dSigmadEprod=0;
119  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
120  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
121 
122 
123  if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile
124 
125  /*G4double Tmax=kinEnergyProj;
126  if (second_part_of_same_type) Tmax = kinEnergyProj/2.;*/
127 
128  G4double E1=kinEnergyProd;
129  G4double E2=kinEnergyProd*1.000001;
130  G4double dE=(E2-E1);
131  G4double sigma1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
132  G4double sigma2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
133 
134  dSigmadEprod=(sigma1-sigma2)/dE;
135  }
136  return dSigmadEprod;
137 
138 
139 
140 }
141 //The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
143 //
144 G4double G4VEmAdjointModel::DiffCrossSectionPerAtomPrimToScatPrim(
145  G4double kinEnergyProj,
146  G4double kinEnergyScatProj,
147  G4double Z,
148  G4double A)
149 { G4double kinEnergyProd = kinEnergyProj - kinEnergyScatProj;
150  G4double dSigmadEprod;
151  if (kinEnergyProd <=0) dSigmadEprod=0;
152  else dSigmadEprod=DiffCrossSectionPerAtomPrimToSecond(kinEnergyProj,kinEnergyProd,Z,A);
153  return dSigmadEprod;
154 
155 }
156 
158 //
159 //The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
160 G4double G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToSecond(
161  const G4Material* aMaterial,
162  G4double kinEnergyProj,
163  G4double kinEnergyProd)
164 {
165  G4double dSigmadEprod=0;
166  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
167  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
168 
169 
170  if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
171  /*G4double Tmax=kinEnergyProj;
172  if (second_part_of_same_type) Tmax = kinEnergyProj/2.;*/
173  G4double E1=kinEnergyProd;
174  G4double E2=kinEnergyProd*1.0001;
175  G4double dE=(E2-E1);
176  G4double sigma1=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E1,1.e20);
177  G4double sigma2=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E2,1.e20);
178  dSigmadEprod=(sigma1-sigma2)/dE;
179  }
180  return dSigmadEprod;
181 
182 
183 
184 }
185 //The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
187 //
188 G4double G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToScatPrim(
189  const G4Material* aMaterial,
190  G4double kinEnergyProj,
191  G4double kinEnergyScatProj)
192 { G4double kinEnergyProd = kinEnergyProj - kinEnergyScatProj;
193  G4double dSigmadEprod;
194  if (kinEnergyProd <=0) dSigmadEprod=0;
195  else dSigmadEprod=DiffCrossSectionPerVolumePrimToSecond(aMaterial,kinEnergyProj,kinEnergyProd);
196  return dSigmadEprod;
197 
198 }
200 //
201 G4double G4VEmAdjointModel::DiffCrossSectionFunction1(G4double kinEnergyProj){
202 
203 
204  G4double bias_factor = CS_biasing_factor*kinEnergyProdForIntegration/kinEnergyProj;
205 
206 
207  if (UseMatrixPerElement ) {
208  return DiffCrossSectionPerAtomPrimToSecond(kinEnergyProj,kinEnergyProdForIntegration,ZSelectedNucleus,ASelectedNucleus)*bias_factor;
209  }
210  else {
211  return DiffCrossSectionPerVolumePrimToSecond(SelectedMaterial,kinEnergyProj,kinEnergyProdForIntegration)*bias_factor;
212  }
213 }
214 
216 //
217 G4double G4VEmAdjointModel::DiffCrossSectionFunction2(G4double kinEnergyProj){
218 
219  G4double bias_factor = CS_biasing_factor*kinEnergyScatProjForIntegration/kinEnergyProj;
220  if (UseMatrixPerElement ) {
221  return DiffCrossSectionPerAtomPrimToScatPrim(kinEnergyProj,kinEnergyScatProjForIntegration,ZSelectedNucleus,ASelectedNucleus)*bias_factor;
222  }
223  else {
224  return DiffCrossSectionPerVolumePrimToScatPrim(SelectedMaterial,kinEnergyProj,kinEnergyScatProjForIntegration)*bias_factor;
225 
226  }
227 
228 }
230 //
231 
232 G4double G4VEmAdjointModel::DiffCrossSectionPerVolumeFunctionForIntegrationOverEkinProj(G4double kinEnergyProd)
233 {
234  return DiffCrossSectionPerVolumePrimToSecond(SelectedMaterial,kinEnergyProjForIntegration,kinEnergyProd);
235 }
237 //
238 std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForSecond(
239  G4double kinEnergyProd,
240  G4double Z,
241  G4double A ,
242  G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
243 {
244  G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
245  ASelectedNucleus= int(A);
246  ZSelectedNucleus=int(Z);
247  kinEnergyProdForIntegration = kinEnergyProd;
248 
249  //compute the vector of integrated cross sections
250  //-------------------
251 
252  G4double minEProj= GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
253  G4double maxEProj= GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
254  G4double E1=minEProj;
255  std::vector< double>* log_ESec_vector = new std::vector< double>();
256  std::vector< double>* log_Prob_vector = new std::vector< double>();
257  log_ESec_vector->clear();
258  log_Prob_vector->clear();
259  log_ESec_vector->push_back(std::log(E1));
260  log_Prob_vector->push_back(-50.);
261 
262  G4double E2=std::pow(10.,double( int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
263  G4double fE=std::pow(10.,1./nbin_pro_decade);
264  G4double int_cross_section=0.;
265 
266  if (std::pow(fE,5.)>(maxEProj/minEProj)) fE = std::pow(maxEProj/minEProj,0.2);
267 
268  while (E1 <maxEProj*0.9999999){
269  //G4cout<<E1<<'\t'<<E2<<G4endl;
270 
271  int_cross_section +=integral.Simpson(this,
272  &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 5);
273  log_ESec_vector->push_back(std::log(std::min(E2,maxEProj)));
274  log_Prob_vector->push_back(std::log(int_cross_section));
275  E1=E2;
276  E2*=fE;
277 
278  }
279  std::vector< std::vector<G4double>* > res_mat;
280  res_mat.clear();
281  if (int_cross_section >0.) {
282  res_mat.push_back(log_ESec_vector);
283  res_mat.push_back(log_Prob_vector);
284  }
285 
286  return res_mat;
287 }
288 
290 //
291 std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForScatProj(
292  G4double kinEnergyScatProj,
293  G4double Z,
294  G4double A ,
295  G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
296 { G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
297  ASelectedNucleus=int(A);
298  ZSelectedNucleus=int(Z);
299  kinEnergyScatProjForIntegration = kinEnergyScatProj;
300 
301  //compute the vector of integrated cross sections
302  //-------------------
303 
304  G4double minEProj= GetSecondAdjEnergyMinForScatProjToProjCase(kinEnergyScatProj);
305  G4double maxEProj= GetSecondAdjEnergyMaxForScatProjToProjCase(kinEnergyScatProj);
306  G4double dEmax=maxEProj-kinEnergyScatProj;
307  G4double dEmin=GetLowEnergyLimit();
308  G4double dE1=dEmin;
309  G4double dE2=dEmin;
310 
311 
312  std::vector< double>* log_ESec_vector = new std::vector< double>();
313  std::vector< double>* log_Prob_vector = new std::vector< double>();
314  log_ESec_vector->push_back(std::log(dEmin));
315  log_Prob_vector->push_back(-50.);
316  G4int nbins=std::max( int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
317  G4double fE=std::pow(dEmax/dEmin,1./nbins);
318 
319 
320 
321 
322 
323  G4double int_cross_section=0.;
324 
325  while (dE1 <dEmax*0.9999999999999){
326  dE2=dE1*fE;
327  int_cross_section +=integral.Simpson(this,
328  &G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 5);
329  //G4cout<<"int_cross_section "<<minEProj+dE1<<'\t'<<int_cross_section<<G4endl;
330  log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj-minEProj)));
331  log_Prob_vector->push_back(std::log(int_cross_section));
332  dE1=dE2;
333 
334  }
335 
336 
337  std::vector< std::vector<G4double> *> res_mat;
338  res_mat.clear();
339  if (int_cross_section >0.) {
340  res_mat.push_back(log_ESec_vector);
341  res_mat.push_back(log_Prob_vector);
342  }
343 
344  return res_mat;
345 }
347 //
348 std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForSecond(
349  G4Material* aMaterial,
350  G4double kinEnergyProd,
351  G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
352 { G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
353  SelectedMaterial= aMaterial;
354  kinEnergyProdForIntegration = kinEnergyProd;
355  //compute the vector of integrated cross sections
356  //-------------------
357 
358  G4double minEProj= GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
359  G4double maxEProj= GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
360  G4double E1=minEProj;
361  std::vector< double>* log_ESec_vector = new std::vector< double>();
362  std::vector< double>* log_Prob_vector = new std::vector< double>();
363  log_ESec_vector->clear();
364  log_Prob_vector->clear();
365  log_ESec_vector->push_back(std::log(E1));
366  log_Prob_vector->push_back(-50.);
367 
368  G4double E2=std::pow(10.,double( int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
369  G4double fE=std::pow(10.,1./nbin_pro_decade);
370  G4double int_cross_section=0.;
371 
372  if (std::pow(fE,5.)>(maxEProj/minEProj)) fE = std::pow(maxEProj/minEProj,0.2);
373 
374  while (E1 <maxEProj*0.9999999){
375 
376  int_cross_section +=integral.Simpson(this,
377  &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 5);
378  log_ESec_vector->push_back(std::log(std::min(E2,maxEProj)));
379  log_Prob_vector->push_back(std::log(int_cross_section));
380  E1=E2;
381  E2*=fE;
382 
383  }
384  std::vector< std::vector<G4double>* > res_mat;
385  res_mat.clear();
386 
387  if (int_cross_section >0.) {
388  res_mat.push_back(log_ESec_vector);
389  res_mat.push_back(log_Prob_vector);
390  }
391 
392 
393 
394  return res_mat;
395 }
396 
398 //
399 std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForScatProj(
400  G4Material* aMaterial,
401  G4double kinEnergyScatProj,
402  G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
403 { G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
404  SelectedMaterial= aMaterial;
405  kinEnergyScatProjForIntegration = kinEnergyScatProj;
406 
407  //compute the vector of integrated cross sections
408  //-------------------
409 
410  G4double minEProj= GetSecondAdjEnergyMinForScatProjToProjCase(kinEnergyScatProj);
411  G4double maxEProj= GetSecondAdjEnergyMaxForScatProjToProjCase(kinEnergyScatProj);
412 
413 
414  G4double dEmax=maxEProj-kinEnergyScatProj;
415  G4double dEmin=GetLowEnergyLimit();
416  G4double dE1=dEmin;
417  G4double dE2=dEmin;
418 
419 
420  std::vector< double>* log_ESec_vector = new std::vector< double>();
421  std::vector< double>* log_Prob_vector = new std::vector< double>();
422  log_ESec_vector->push_back(std::log(dEmin));
423  log_Prob_vector->push_back(-50.);
424  G4int nbins=std::max( int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
425  G4double fE=std::pow(dEmax/dEmin,1./nbins);
426 
427  G4double int_cross_section=0.;
428 
429  while (dE1 <dEmax*0.9999999999999){
430  dE2=dE1*fE;
431  int_cross_section +=integral.Simpson(this,
432  &G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 5);
433  log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj-minEProj)));
434  log_Prob_vector->push_back(std::log(int_cross_section));
435  dE1=dE2;
436 
437  }
438 
439 
440 
441 
442 
443  std::vector< std::vector<G4double> *> res_mat;
444  res_mat.clear();
445  if (int_cross_section >0.) {
446  res_mat.push_back(log_ESec_vector);
447  res_mat.push_back(log_Prob_vector);
448  }
449 
450  return res_mat;
451 }
453 //
454 G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex,G4double aPrimEnergy,G4bool IsScatProjToProjCase)
455 {
456 
457 
458  G4AdjointCSMatrix* theMatrix= (*pOnCSMatrixForProdToProjBackwardScattering)[MatrixIndex];
459  if (IsScatProjToProjCase) theMatrix= (*pOnCSMatrixForScatProjToProjBackwardScattering)[MatrixIndex];
460  std::vector< double>* theLogPrimEnergyVector = theMatrix->GetLogPrimEnergyVector();
461 
462  if (theLogPrimEnergyVector->size() ==0){
463  G4cout<<"No data are contained in the given AdjointCSMatrix!"<<G4endl;
464  G4cout<<"The sampling procedure will be stopped."<<G4endl;
465  return 0.;
466 
467  }
468 
469  G4AdjointInterpolator* theInterpolator=G4AdjointInterpolator::GetInstance();
470  G4double aLogPrimEnergy = std::log(aPrimEnergy);
471  size_t ind =theInterpolator->FindPositionForLogVector(aLogPrimEnergy,*theLogPrimEnergyVector);
472 
473 
474  G4double aLogPrimEnergy1,aLogPrimEnergy2;
475  G4double aLogCS1,aLogCS2;
476  G4double log01,log02;
477  std::vector< double>* aLogSecondEnergyVector1 =0;
478  std::vector< double>* aLogSecondEnergyVector2 =0;
479  std::vector< double>* aLogProbVector1=0;
480  std::vector< double>* aLogProbVector2=0;
481  std::vector< size_t>* aLogProbVectorIndex1=0;
482  std::vector< size_t>* aLogProbVectorIndex2=0;
483 
484  theMatrix->GetData(ind, aLogPrimEnergy1,aLogCS1,log01, aLogSecondEnergyVector1,aLogProbVector1,aLogProbVectorIndex1);
485  theMatrix->GetData(ind+1, aLogPrimEnergy2,aLogCS2,log02, aLogSecondEnergyVector2,aLogProbVector2,aLogProbVectorIndex2);
486 
487  G4double rand_var = G4UniformRand();
488  G4double log_rand_var= std::log(rand_var);
489  G4double log_Tcut =std::log(currentTcutForDirectSecond);
490  G4double Esec=0;
491  G4double log_dE1,log_dE2;
492  G4double log_rand_var1,log_rand_var2;
493  G4double log_E1,log_E2;
494  log_rand_var1=log_rand_var;
495  log_rand_var2=log_rand_var;
496 
497  G4double Emin=0.;
498  G4double Emax=0.;
499  if (theMatrix->IsScatProjToProjCase()){ //case where Tcut plays a role
500  Emin=GetSecondAdjEnergyMinForScatProjToProjCase(aPrimEnergy,currentTcutForDirectSecond);
501  Emax=GetSecondAdjEnergyMaxForScatProjToProjCase(aPrimEnergy);
502  G4double dE=0;
503  if (Emin < Emax ){
504  if (ApplyCutInRange) {
505  if (second_part_of_same_type && currentTcutForDirectSecond>aPrimEnergy) return aPrimEnergy;
506 
507  log_rand_var1=log_rand_var+theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector1,*aLogProbVector1);
508  log_rand_var2=log_rand_var+theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector2,*aLogProbVector2);
509 
510  }
511  log_dE1 = theInterpolator->Interpolate(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
512  log_dE2 = theInterpolator->Interpolate(log_rand_var2,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
513  dE=std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_dE1,log_dE2));
514  }
515 
516  Esec = aPrimEnergy +dE;
517  Esec=std::max(Esec,Emin);
518  Esec=std::min(Esec,Emax);
519 
520  }
521  else { //Tcut condition is already full-filled
522 
523  log_E1 = theInterpolator->Interpolate(log_rand_var,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
524  log_E2 = theInterpolator->Interpolate(log_rand_var,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
525 
526  Esec = std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_E1,log_E2));
527  Emin=GetSecondAdjEnergyMinForProdToProjCase(aPrimEnergy);
528  Emax=GetSecondAdjEnergyMaxForProdToProjCase(aPrimEnergy);
529  Esec=std::max(Esec,Emin);
530  Esec=std::min(Esec,Emax);
531 
532  }
533 
534  return Esec;
535 
536 
537 
538 
539 
540 }
541 
543 //
544 G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(G4double aPrimEnergy,G4bool IsScatProjToProjCase)
545 { SelectCSMatrix(IsScatProjToProjCase);
546  return SampleAdjSecEnergyFromCSMatrix(indexOfUsedCrossSectionMatrix, aPrimEnergy, IsScatProjToProjCase);
547 }
549 //
550 void G4VEmAdjointModel::SelectCSMatrix(G4bool IsScatProjToProjCase)
551 {
552  indexOfUsedCrossSectionMatrix=0;
553  if (!UseMatrixPerElement) indexOfUsedCrossSectionMatrix = currentMaterialIndex;
554  else if (!UseOnlyOneMatrixForAllElements) { //Select Material
555  std::vector<G4double>* CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
556  lastCS=lastAdjointCSForScatProjToProjCase;
557  if ( !IsScatProjToProjCase) {
558  CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
559  lastCS=lastAdjointCSForProdToProjCase;
560  }
561  G4double rand_var= G4UniformRand();
562  G4double SumCS=0.;
563  size_t ind=0;
564  for (size_t i=0;i<CS_Vs_Element->size();i++){
565  SumCS+=(*CS_Vs_Element)[i];
566  if (rand_var<=SumCS/lastCS){
567  ind=i;
568  break;
569  }
570  }
571  indexOfUsedCrossSectionMatrix = currentMaterial->GetElement(ind)->GetIndex();
572  }
573 }
575 //
576 G4double G4VEmAdjointModel::SampleAdjSecEnergyFromDiffCrossSectionPerAtom(G4double prim_energy,G4bool IsScatProjToProjCase)
577 {
578  // here we try to use the rejection method
579  //-----------------------------------------
580 
581  G4double E=0;
582  G4double x,xmin,greject,q;
583  if ( IsScatProjToProjCase){
584  G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(prim_energy);
585  G4double Emin= prim_energy+currentTcutForDirectSecond;
586  xmin=Emin/Emax;
587  G4double grejmax = DiffCrossSectionPerAtomPrimToScatPrim(Emin,prim_energy,1)*prim_energy;
588 
589  do {
590  q = G4UniformRand();
591  x = 1./(q*(1./xmin -1.) +1.);
592  E=x*Emax;
593  greject = DiffCrossSectionPerAtomPrimToScatPrim( E,prim_energy ,1)*prim_energy;
594 
595  }
596 
597  while( greject < G4UniformRand()*grejmax );
598 
599  }
600  else {
601  G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(prim_energy);
602  G4double Emin= GetSecondAdjEnergyMinForProdToProjCase(prim_energy);;
603  xmin=Emin/Emax;
604  G4double grejmax = DiffCrossSectionPerAtomPrimToSecond(Emin,prim_energy,1);
605  do {
606  q = G4UniformRand();
607  x = std::pow(xmin, q);
608  E=x*Emax;
609  greject = DiffCrossSectionPerAtomPrimToSecond( E,prim_energy ,1);
610 
611  }
612 
613  while( greject < G4UniformRand()*grejmax );
614 
615 
616 
617  }
618 
619  return E;
620 }
621 
623 //
624 void G4VEmAdjointModel::CorrectPostStepWeight(G4ParticleChange* fParticleChange,
625  G4double old_weight,
626  G4double adjointPrimKinEnergy,
627  G4double projectileKinEnergy,
628  G4bool IsScatProjToProjCase)
629 {
630  G4double new_weight=old_weight;
631  G4double w_corr =1./CS_biasing_factor;
632  w_corr*=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
633 
634 
635  lastCS=lastAdjointCSForScatProjToProjCase;
636  if ( !IsScatProjToProjCase) lastCS=lastAdjointCSForProdToProjCase;
637  if ((adjointPrimKinEnergy-preStepEnergy)/preStepEnergy>0.001){ //Is that in all cases needed???
638  G4double post_stepCS=AdjointCrossSection(currentCouple, adjointPrimKinEnergy
639  ,IsScatProjToProjCase );
640  if (post_stepCS>0 && lastCS>0) w_corr*=post_stepCS/lastCS;
641  }
642 
643  new_weight*=w_corr;
644 
645  //G4cout<<"Post step "<<new_weight<<'\t'<<w_corr<<'\t'<<old_weight<<G4endl;
646  new_weight*=projectileKinEnergy/adjointPrimKinEnergy;//This is needed due to the biasing of diff CS
647  //by the factor adjointPrimKinEnergy/projectileKinEnergy
648 
649 
650 
651  fParticleChange->SetParentWeightByProcess(false);
652  fParticleChange->SetSecondaryWeightByProcess(false);
653  fParticleChange->ProposeParentWeight(new_weight);
654 }
656 //
657 G4double G4VEmAdjointModel::GetSecondAdjEnergyMaxForScatProjToProjCase(G4double kinEnergyScatProj)
658 { G4double maxEProj= HighEnergyLimit;
659  if (second_part_of_same_type) maxEProj=std::min(kinEnergyScatProj*2.,HighEnergyLimit);
660  return maxEProj;
661 }
663 //
664 G4double G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy,G4double Tcut)
665 { G4double Emin=PrimAdjEnergy;
666  if (ApplyCutInRange) Emin=PrimAdjEnergy+Tcut;
667  return Emin;
668 }
670 //
671 G4double G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProjCase(G4double )
672 { return HighEnergyLimit;
673 }
675 //
676 G4double G4VEmAdjointModel::GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
677 { G4double minEProj=PrimAdjEnergy;
678  if (second_part_of_same_type) minEProj=PrimAdjEnergy*2.;
679  return minEProj;
680 }
682 //
683 void G4VEmAdjointModel::DefineCurrentMaterial(const G4MaterialCutsCouple* couple)
684 { if(couple != currentCouple) {
685  currentCouple = const_cast<G4MaterialCutsCouple*> (couple);
686  currentMaterial = const_cast<G4Material*> (couple->GetMaterial());
687  currentCoupleIndex = couple->GetIndex();
688  currentMaterialIndex = currentMaterial->GetIndex();
689  size_t idx=56;
690  currentTcutForDirectSecond =0.00000000001;
691  if (theAdjEquivOfDirectSecondPartDef) {
692  if (theAdjEquivOfDirectSecondPartDef == G4AdjointGamma::AdjointGamma()) idx = 0;
693  else if (theAdjEquivOfDirectSecondPartDef == G4AdjointElectron::AdjointElectron()) idx = 1;
694  else if (theAdjEquivOfDirectSecondPartDef == G4AdjointPositron::AdjointPositron()) idx = 2;
695  if (idx <56){
696  const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
697  currentTcutForDirectSecond=(*aVec)[currentCoupleIndex];
698  }
699  }
700 
701 
702  }
703 }
705 //
706 void G4VEmAdjointModel::SetHighEnergyLimit(G4double aVal)
707 { HighEnergyLimit=aVal;
708  if (theDirectEMModel) theDirectEMModel->SetHighEnergyLimit( aVal);
709 }
711 //
712 void G4VEmAdjointModel::SetLowEnergyLimit(G4double aVal)
713 {
714  LowEnergyLimit=aVal;
715  if (theDirectEMModel) theDirectEMModel->SetLowEnergyLimit( aVal);
716 }
718 //
719 void G4VEmAdjointModel::SetAdjointEquivalentOfDirectPrimaryParticleDefinition(G4ParticleDefinition* aPart)
720 {
721  theAdjEquivOfDirectPrimPartDef=aPart;
722  if (theAdjEquivOfDirectPrimPartDef->GetParticleName() =="adj_e-")
723  theDirectPrimaryPartDef=G4Electron::Electron();
724  if (theAdjEquivOfDirectPrimPartDef->GetParticleName() =="adj_gamma")
725  theDirectPrimaryPartDef=G4Gamma::Gamma();
726 }
Generated on Sat May 25 2013 14:33:29 for Geant4 by   doxygen 1.8.4