ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4VEmAdjointModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4VEmAdjointModel.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 "G4VEmAdjointModel.hh"
28 #include "G4AdjointCSManager.hh"
29 #include "G4Integrator.hh"
30 #include "G4TrackStatus.hh"
31 #include "G4ParticleChange.hh"
32 #include "G4AdjointElectron.hh"
33 #include "G4AdjointGamma.hh"
34 #include "G4AdjointPositron.hh"
35 #include "G4AdjointInterpolator.hh"
36 #include "G4PhysicsTable.hh"
37 
39 //
40 G4VEmAdjointModel::G4VEmAdjointModel(const G4String& nam):
41 name(nam)
42 // lowLimit(0.1*keV), highLimit(100.0*TeV), fluc(0), name(nam), pParticleChange(0)
43 {
44  model_index = G4AdjointCSManager::GetAdjointCSManager()->RegisterEmAdjointModel(this);
45  second_part_of_same_type =false;
46  theDirectEMModel=0;
47  mass_ratio_product=1.;
48  mass_ratio_projectile=1.;
49  currentCouple=0;
50  additional_weight_correction_factor_for_post_step_outside_model=1.;
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  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
269  while (E1 <maxEProj*0.9999999){
270  //G4cout<<E1<<'\t'<<E2<<G4endl;
271 
272  int_cross_section +=integral.Simpson(this,
273  &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 5);
274  log_ESec_vector->push_back(std::log(std::min(E2,maxEProj)));
275  log_Prob_vector->push_back(std::log(int_cross_section));
276  E1=E2;
277  E2*=fE;
278 
279  }
280  std::vector< std::vector<G4double>* > res_mat;
281  res_mat.clear();
282  if (int_cross_section >0.) {
283  res_mat.push_back(log_ESec_vector);
284  res_mat.push_back(log_Prob_vector);
285  }
286 
287  return res_mat;
288 }
289 
291 //
292 std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForScatProj(
293  G4double kinEnergyScatProj,
294  G4double Z,
295  G4double A ,
296  G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
297 { G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
298  ASelectedNucleus=int(A);
299  ZSelectedNucleus=int(Z);
300  kinEnergyScatProjForIntegration = kinEnergyScatProj;
301 
302  //compute the vector of integrated cross sections
303  //-------------------
304 
305  G4double minEProj= GetSecondAdjEnergyMinForScatProjToProjCase(kinEnergyScatProj);
306  G4double maxEProj= GetSecondAdjEnergyMaxForScatProjToProjCase(kinEnergyScatProj);
307  G4double dEmax=maxEProj-kinEnergyScatProj;
308  G4double dEmin=GetLowEnergyLimit();
309  G4double dE1=dEmin;
310  G4double dE2=dEmin;
311 
312 
313  std::vector< double>* log_ESec_vector = new std::vector< double>();
314  std::vector< double>* log_Prob_vector = new std::vector< double>();
315  log_ESec_vector->push_back(std::log(dEmin));
316  log_Prob_vector->push_back(-50.);
317  G4int nbins=std::max( int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
318  G4double fE=std::pow(dEmax/dEmin,1./nbins);
319 
320 
321 
322 
323 
324  G4double int_cross_section=0.;
325 
326  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
327  while (dE1 <dEmax*0.9999999999999){
328  dE2=dE1*fE;
329  int_cross_section +=integral.Simpson(this,
330  &G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 5);
331  //G4cout<<"int_cross_section "<<minEProj+dE1<<'\t'<<int_cross_section<<G4endl;
332  log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj-minEProj)));
333  log_Prob_vector->push_back(std::log(int_cross_section));
334  dE1=dE2;
335 
336  }
337 
338 
339  std::vector< std::vector<G4double> *> res_mat;
340  res_mat.clear();
341  if (int_cross_section >0.) {
342  res_mat.push_back(log_ESec_vector);
343  res_mat.push_back(log_Prob_vector);
344  }
345 
346  return res_mat;
347 }
349 //
350 std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForSecond(
351  G4Material* aMaterial,
352  G4double kinEnergyProd,
353  G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
354 { G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
355  SelectedMaterial= aMaterial;
356  kinEnergyProdForIntegration = kinEnergyProd;
357  //compute the vector of integrated cross sections
358  //-------------------
359 
360  G4double minEProj= GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
361  G4double maxEProj= GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
362  G4double E1=minEProj;
363  std::vector< double>* log_ESec_vector = new std::vector< double>();
364  std::vector< double>* log_Prob_vector = new std::vector< double>();
365  log_ESec_vector->clear();
366  log_Prob_vector->clear();
367  log_ESec_vector->push_back(std::log(E1));
368  log_Prob_vector->push_back(-50.);
369 
370  G4double E2=std::pow(10.,double( int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
371  G4double fE=std::pow(10.,1./nbin_pro_decade);
372  G4double int_cross_section=0.;
373 
374  if (std::pow(fE,5.)>(maxEProj/minEProj)) fE = std::pow(maxEProj/minEProj,0.2);
375 
376  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
377  while (E1 <maxEProj*0.9999999){
378 
379  int_cross_section +=integral.Simpson(this,
380  &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 5);
381  log_ESec_vector->push_back(std::log(std::min(E2,maxEProj)));
382  log_Prob_vector->push_back(std::log(int_cross_section));
383  E1=E2;
384  E2*=fE;
385 
386  }
387  std::vector< std::vector<G4double>* > res_mat;
388  res_mat.clear();
389 
390  if (int_cross_section >0.) {
391  res_mat.push_back(log_ESec_vector);
392  res_mat.push_back(log_Prob_vector);
393  }
394 
395 
396 
397  return res_mat;
398 }
399 
401 //
402 std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForScatProj(
403  G4Material* aMaterial,
404  G4double kinEnergyScatProj,
405  G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
406 { G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
407  SelectedMaterial= aMaterial;
408  kinEnergyScatProjForIntegration = kinEnergyScatProj;
409 
410  //compute the vector of integrated cross sections
411  //-------------------
412 
413  G4double minEProj= GetSecondAdjEnergyMinForScatProjToProjCase(kinEnergyScatProj);
414  G4double maxEProj= GetSecondAdjEnergyMaxForScatProjToProjCase(kinEnergyScatProj);
415 
416 
417  G4double dEmax=maxEProj-kinEnergyScatProj;
418  G4double dEmin=GetLowEnergyLimit();
419  G4double dE1=dEmin;
420  G4double dE2=dEmin;
421 
422 
423  std::vector< double>* log_ESec_vector = new std::vector< double>();
424  std::vector< double>* log_Prob_vector = new std::vector< double>();
425  log_ESec_vector->push_back(std::log(dEmin));
426  log_Prob_vector->push_back(-50.);
427  G4int nbins=std::max( int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
428  G4double fE=std::pow(dEmax/dEmin,1./nbins);
429 
430  G4double int_cross_section=0.;
431 
432  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
433  while (dE1 <dEmax*0.9999999999999){
434  dE2=dE1*fE;
435  int_cross_section +=integral.Simpson(this,
436  &G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 5);
437  log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj-minEProj)));
438  log_Prob_vector->push_back(std::log(int_cross_section));
439  dE1=dE2;
440 
441  }
442 
443 
444 
445 
446 
447  std::vector< std::vector<G4double> *> res_mat;
448  res_mat.clear();
449  if (int_cross_section >0.) {
450  res_mat.push_back(log_ESec_vector);
451  res_mat.push_back(log_Prob_vector);
452  }
453 
454  return res_mat;
455 }
457 //
458 G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex,G4double aPrimEnergy,G4bool IsScatProjToProjCase)
459 {
460 
461 
462  G4AdjointCSMatrix* theMatrix= (*pOnCSMatrixForProdToProjBackwardScattering)[MatrixIndex];
463  if (IsScatProjToProjCase) theMatrix= (*pOnCSMatrixForScatProjToProjBackwardScattering)[MatrixIndex];
464  std::vector< double>* theLogPrimEnergyVector = theMatrix->GetLogPrimEnergyVector();
465 
466  if (theLogPrimEnergyVector->size() ==0){
467  G4cout<<"No data are contained in the given AdjointCSMatrix!"<<G4endl;
468  G4cout<<"The sampling procedure will be stopped."<<G4endl;
469  return 0.;
470 
471  }
472 
473  G4AdjointInterpolator* theInterpolator=G4AdjointInterpolator::GetInstance();
474  G4double aLogPrimEnergy = std::log(aPrimEnergy);
475  size_t ind =theInterpolator->FindPositionForLogVector(aLogPrimEnergy,*theLogPrimEnergyVector);
476 
477 
478  G4double aLogPrimEnergy1,aLogPrimEnergy2;
479  G4double aLogCS1,aLogCS2;
480  G4double log01,log02;
481  std::vector< double>* aLogSecondEnergyVector1 =0;
482  std::vector< double>* aLogSecondEnergyVector2 =0;
483  std::vector< double>* aLogProbVector1=0;
484  std::vector< double>* aLogProbVector2=0;
485  std::vector< size_t>* aLogProbVectorIndex1=0;
486  std::vector< size_t>* aLogProbVectorIndex2=0;
487 
488  theMatrix->GetData(ind, aLogPrimEnergy1,aLogCS1,log01, aLogSecondEnergyVector1,aLogProbVector1,aLogProbVectorIndex1);
489  theMatrix->GetData(ind+1, aLogPrimEnergy2,aLogCS2,log02, aLogSecondEnergyVector2,aLogProbVector2,aLogProbVectorIndex2);
490 
491  G4double rand_var = G4UniformRand();
492  G4double log_rand_var= std::log(rand_var);
493  G4double log_Tcut =std::log(currentTcutForDirectSecond);
494  G4double Esec=0;
495  G4double log_dE1,log_dE2;
496  G4double log_rand_var1,log_rand_var2;
497  G4double log_E1,log_E2;
498  log_rand_var1=log_rand_var;
499  log_rand_var2=log_rand_var;
500 
501  G4double Emin=0.;
502  G4double Emax=0.;
503  if (theMatrix->IsScatProjToProjCase()){ //case where Tcut plays a role
504  Emin=GetSecondAdjEnergyMinForScatProjToProjCase(aPrimEnergy,currentTcutForDirectSecond);
505  Emax=GetSecondAdjEnergyMaxForScatProjToProjCase(aPrimEnergy);
506  G4double dE=0;
507  if (Emin < Emax ){
508  if (ApplyCutInRange) {
509  if (second_part_of_same_type && currentTcutForDirectSecond>aPrimEnergy) return aPrimEnergy;
510 
511  log_rand_var1=log_rand_var+theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector1,*aLogProbVector1);
512  log_rand_var2=log_rand_var+theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector2,*aLogProbVector2);
513 
514  }
515  log_dE1 = theInterpolator->Interpolate(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
516  log_dE2 = theInterpolator->Interpolate(log_rand_var2,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
517  dE=std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_dE1,log_dE2));
518  }
519 
520  Esec = aPrimEnergy +dE;
521  Esec=std::max(Esec,Emin);
522  Esec=std::min(Esec,Emax);
523 
524  }
525  else { //Tcut condition is already full-filled
526 
527  log_E1 = theInterpolator->Interpolate(log_rand_var,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
528  log_E2 = theInterpolator->Interpolate(log_rand_var,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
529 
530  Esec = std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_E1,log_E2));
531  Emin=GetSecondAdjEnergyMinForProdToProjCase(aPrimEnergy);
532  Emax=GetSecondAdjEnergyMaxForProdToProjCase(aPrimEnergy);
533  Esec=std::max(Esec,Emin);
534  Esec=std::min(Esec,Emax);
535 
536  }
537 
538  return Esec;
539 
540 
541 
542 
543 
544 }
545 
547 //
548 G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(G4double aPrimEnergy,G4bool IsScatProjToProjCase)
549 { SelectCSMatrix(IsScatProjToProjCase);
550  return SampleAdjSecEnergyFromCSMatrix(indexOfUsedCrossSectionMatrix, aPrimEnergy, IsScatProjToProjCase);
551 }
553 //
554 void G4VEmAdjointModel::SelectCSMatrix(G4bool IsScatProjToProjCase)
555 {
556  indexOfUsedCrossSectionMatrix=0;
557  if (!UseMatrixPerElement) indexOfUsedCrossSectionMatrix = currentMaterialIndex;
558  else if (!UseOnlyOneMatrixForAllElements) { //Select Material
559  std::vector<G4double>* CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
560  lastCS=lastAdjointCSForScatProjToProjCase;
561  if ( !IsScatProjToProjCase) {
562  CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
563  lastCS=lastAdjointCSForProdToProjCase;
564  }
565  G4double rand_var= G4UniformRand();
566  G4double SumCS=0.;
567  size_t ind=0;
568  for (size_t i=0;i<CS_Vs_Element->size();i++){
569  SumCS+=(*CS_Vs_Element)[i];
570  if (rand_var<=SumCS/lastCS){
571  ind=i;
572  break;
573  }
574  }
575  indexOfUsedCrossSectionMatrix = currentMaterial->GetElement(ind)->GetIndex();
576  }
577 }
579 //
580 G4double G4VEmAdjointModel::SampleAdjSecEnergyFromDiffCrossSectionPerAtom(G4double prim_energy,G4bool IsScatProjToProjCase)
581 {
582  // here we try to use the rejection method
583  //-----------------------------------------
584 
585  const G4int iimax = 1000;
586  G4double E=0;
587  G4double x,xmin,greject,q;
588  if ( IsScatProjToProjCase){
589  G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(prim_energy);
590  G4double Emin= prim_energy+currentTcutForDirectSecond;
591  xmin=Emin/Emax;
592  G4double grejmax = DiffCrossSectionPerAtomPrimToScatPrim(Emin,prim_energy,1)*prim_energy;
593 
594  G4int ii =0;
595  do {
596  q = G4UniformRand();
597  x = 1./(q*(1./xmin -1.) +1.);
598  E=x*Emax;
599  greject = DiffCrossSectionPerAtomPrimToScatPrim( E,prim_energy ,1)*prim_energy;
600  ++ii;
601  if(ii >= iimax) { break; }
602  }
603  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
604  while( greject < G4UniformRand()*grejmax );
605 
606  }
607  else {
608  G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(prim_energy);
609  G4double Emin= GetSecondAdjEnergyMinForProdToProjCase(prim_energy);;
610  xmin=Emin/Emax;
611  G4double grejmax = DiffCrossSectionPerAtomPrimToSecond(Emin,prim_energy,1);
612  G4int ii =0;
613  do {
614  q = G4UniformRand();
615  x = std::pow(xmin, q);
616  E=x*Emax;
617  greject = DiffCrossSectionPerAtomPrimToSecond( E,prim_energy ,1);
618  ++ii;
619  if(ii >= iimax) { break; }
620  }
621  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
622  while( greject < G4UniformRand()*grejmax );
623 
624 
625 
626  }
627 
628  return E;
629 }
630 
632 //
633 void G4VEmAdjointModel::CorrectPostStepWeight(G4ParticleChange* fParticleChange,
634  G4double old_weight,
635  G4double adjointPrimKinEnergy,
636  G4double projectileKinEnergy,
637  G4bool IsScatProjToProjCase)
638 {
639  G4double new_weight=old_weight;
640  G4double w_corr =1./CS_biasing_factor;
641  w_corr*=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
642 
643 
644  lastCS=lastAdjointCSForScatProjToProjCase;
645  if ( !IsScatProjToProjCase) lastCS=lastAdjointCSForProdToProjCase;
646  if ((adjointPrimKinEnergy-preStepEnergy)/preStepEnergy>0.001){ //Is that in all cases needed???
647  G4double post_stepCS=AdjointCrossSection(currentCouple, adjointPrimKinEnergy
648  ,IsScatProjToProjCase );
649  if (post_stepCS>0 && lastCS>0) w_corr*=post_stepCS/lastCS;
650  }
651 
652  new_weight*=w_corr;
653 
654  //G4cout<<"Post step "<<new_weight<<'\t'<<w_corr<<'\t'<<old_weight<<G4endl;
655  new_weight*=projectileKinEnergy/adjointPrimKinEnergy;//This is needed due to the biasing of diff CS
656  //by the factor adjointPrimKinEnergy/projectileKinEnergy
657 
658 
659 
660  fParticleChange->SetParentWeightByProcess(false);
661  fParticleChange->SetSecondaryWeightByProcess(false);
662  fParticleChange->ProposeParentWeight(new_weight);
663 }
665 //
666 G4double G4VEmAdjointModel::GetSecondAdjEnergyMaxForScatProjToProjCase(G4double kinEnergyScatProj)
667 { G4double maxEProj= HighEnergyLimit;
668  if (second_part_of_same_type) maxEProj=std::min(kinEnergyScatProj*2.,HighEnergyLimit);
669  return maxEProj;
670 }
672 //
673 G4double G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy,G4double Tcut)
674 { G4double Emin=PrimAdjEnergy;
675  if (ApplyCutInRange) Emin=PrimAdjEnergy+Tcut;
676  return Emin;
677 }
679 //
680 G4double G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProjCase(G4double )
681 { return HighEnergyLimit;
682 }
684 //
685 G4double G4VEmAdjointModel::GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
686 { G4double minEProj=PrimAdjEnergy;
687  if (second_part_of_same_type) minEProj=PrimAdjEnergy*2.;
688  return minEProj;
689 }
691 //
692 void G4VEmAdjointModel::DefineCurrentMaterial(const G4MaterialCutsCouple* couple)
693 { if(couple != currentCouple) {
694  currentCouple = const_cast<G4MaterialCutsCouple*> (couple);
695  currentMaterial = const_cast<G4Material*> (couple->GetMaterial());
696  currentCoupleIndex = couple->GetIndex();
697  currentMaterialIndex = currentMaterial->GetIndex();
698  size_t idx=56;
699  currentTcutForDirectSecond =0.00000000001;
700  if (theAdjEquivOfDirectSecondPartDef) {
701  if (theAdjEquivOfDirectSecondPartDef == G4AdjointGamma::AdjointGamma()) idx = 0;
702  else if (theAdjEquivOfDirectSecondPartDef == G4AdjointElectron::AdjointElectron()) idx = 1;
703  else if (theAdjEquivOfDirectSecondPartDef == G4AdjointPositron::AdjointPositron()) idx = 2;
704  if (idx <56){
705  const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
706  currentTcutForDirectSecond=(*aVec)[currentCoupleIndex];
707  }
708  }
709 
710 
711  }
712 }
714 //
715 void G4VEmAdjointModel::SetHighEnergyLimit(G4double aVal)
716 { HighEnergyLimit=aVal;
717  if (theDirectEMModel) theDirectEMModel->SetHighEnergyLimit( aVal);
718 }
720 //
721 void G4VEmAdjointModel::SetLowEnergyLimit(G4double aVal)
722 {
723  LowEnergyLimit=aVal;
724  if (theDirectEMModel) theDirectEMModel->SetLowEnergyLimit( aVal);
725 }
727 //
728 void G4VEmAdjointModel::SetAdjointEquivalentOfDirectPrimaryParticleDefinition(G4ParticleDefinition* aPart)
729 {
730  theAdjEquivOfDirectPrimPartDef=aPart;
731  if (theAdjEquivOfDirectPrimPartDef->GetParticleName() =="adj_e-")
732  theDirectPrimaryPartDef=G4Electron::Electron();
733  if (theAdjEquivOfDirectPrimPartDef->GetParticleName() =="adj_gamma")
734  theDirectPrimaryPartDef=G4Gamma::Gamma();
735 }