Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
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 // $Id: G4VEmAdjointModel.cc 100341 2016-10-18 08:02:25Z 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 //
42 name(nam)
43 // lowLimit(0.1*keV), highLimit(100.0*TeV), fluc(0), name(nam), pParticleChange(0)
44 {
50  currentCouple=0;
52 }
54 //
56 {;}
58 //
60  G4double primEnergy,
61  G4bool IsScatProjToProjCase)
62 {
63  DefineCurrentMaterial(aCouple);
64  preStepEnergy=primEnergy;
65 
66  std::vector<G4double>* CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
67  if (IsScatProjToProjCase) CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
69  this,
70  primEnergy,
72  IsScatProjToProjCase,
73  *CS_Vs_Element);
74  if (IsScatProjToProjCase) lastAdjointCSForScatProjToProjCase = lastCS;
76 
77 
78 
79  return lastCS;
80 
81 }
83 //
85  G4double primEnergy,
86  G4bool IsScatProjToProjCase)
87 {
88  return AdjointCrossSection(aCouple, primEnergy,
89  IsScatProjToProjCase);
90 
91  /*
92  //To continue
93  DefineCurrentMaterial(aCouple);
94  preStepEnergy=primEnergy;
95  if (IsScatProjToProjCase){
96  G4double ekin=primEnergy*mass_ratio_projectile;
97  lastCS = G4AdjointCSManager::GetAdjointCSManager()->GetAdjointSigma(ekin, model_index,true, aCouple);
98  lastAdjointCSForScatProjToProjCase = lastCS;
99  //G4cout<<ekin<<std::endl;
100  }
101  else {
102  G4double ekin=primEnergy*mass_ratio_product;
103  lastCS = G4AdjointCSManager::GetAdjointCSManager()->GetAdjointSigma(ekin, model_index,false, aCouple);
104  lastAdjointCSForProdToProjCase = lastCS;
105  //G4cout<<ekin<<std::endl;
106  }
107  return lastCS;
108  */
109 }
111 //
112 //General implementation correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
114  G4double kinEnergyProj,
115  G4double kinEnergyProd,
116  G4double Z,
117  G4double A)
118 {
119  G4double dSigmadEprod=0;
120  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
121  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
122 
123 
124  if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile
125 
126  /*G4double Tmax=kinEnergyProj;
127  if (second_part_of_same_type) Tmax = kinEnergyProj/2.;*/
128 
129  G4double E1=kinEnergyProd;
130  G4double E2=kinEnergyProd*1.000001;
131  G4double dE=(E2-E1);
134 
135  dSigmadEprod=(sigma1-sigma2)/dE;
136  }
137  return dSigmadEprod;
138 
139 
140 
141 }
142 //The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
144 //
146  G4double kinEnergyProj,
147  G4double kinEnergyScatProj,
148  G4double Z,
149  G4double A)
150 { G4double kinEnergyProd = kinEnergyProj - kinEnergyScatProj;
151  G4double dSigmadEprod;
152  if (kinEnergyProd <=0) dSigmadEprod=0;
153  else dSigmadEprod=DiffCrossSectionPerAtomPrimToSecond(kinEnergyProj,kinEnergyProd,Z,A);
154  return dSigmadEprod;
155 
156 }
157 
159 //
160 //The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
162  const G4Material* aMaterial,
163  G4double kinEnergyProj,
164  G4double kinEnergyProd)
165 {
166  G4double dSigmadEprod=0;
167  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
168  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
169 
170 
171  if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
172  /*G4double Tmax=kinEnergyProj;
173  if (second_part_of_same_type) Tmax = kinEnergyProj/2.;*/
174  G4double E1=kinEnergyProd;
175  G4double E2=kinEnergyProd*1.0001;
176  G4double dE=(E2-E1);
177  G4double sigma1=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E1,1.e20);
178  G4double sigma2=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E2,1.e20);
179  dSigmadEprod=(sigma1-sigma2)/dE;
180  }
181  return dSigmadEprod;
182 
183 
184 
185 }
186 //The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
188 //
190  const G4Material* aMaterial,
191  G4double kinEnergyProj,
192  G4double kinEnergyScatProj)
193 { G4double kinEnergyProd = kinEnergyProj - kinEnergyScatProj;
194  G4double dSigmadEprod;
195  if (kinEnergyProd <=0) dSigmadEprod=0;
196  else dSigmadEprod=DiffCrossSectionPerVolumePrimToSecond(aMaterial,kinEnergyProj,kinEnergyProd);
197  return dSigmadEprod;
198 
199 }
201 //
203 
204 
205  G4double bias_factor = CS_biasing_factor*kinEnergyProdForIntegration/kinEnergyProj;
206 
207 
208  if (UseMatrixPerElement ) {
210  }
211  else {
213  }
214 }
215 
217 //
219 
220  G4double bias_factor = CS_biasing_factor*kinEnergyScatProjForIntegration/kinEnergyProj;
221  if (UseMatrixPerElement ) {
223  }
224  else {
226 
227  }
228 
229 }
231 //
232 
234 {
236 }
238 //
240  G4double kinEnergyProd,
241  G4double Z,
242  G4double A ,
243  G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
244 {
246  ASelectedNucleus= int(A);
248  kinEnergyProdForIntegration = kinEnergyProd;
249 
250  //compute the vector of integrated cross sections
251  //-------------------
252 
253  G4double minEProj= GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
254  G4double maxEProj= GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
255  G4double E1=minEProj;
256  std::vector< double>* log_ESec_vector = new std::vector< double>();
257  std::vector< double>* log_Prob_vector = new std::vector< double>();
258  log_ESec_vector->clear();
259  log_Prob_vector->clear();
260  log_ESec_vector->push_back(std::log(E1));
261  log_Prob_vector->push_back(-50.);
262 
263  G4double E2=std::pow(10.,double( int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
264  G4double fE=std::pow(10.,1./nbin_pro_decade);
265  G4double int_cross_section=0.;
266 
267  if (std::pow(fE,5.)>(maxEProj/minEProj)) fE = std::pow(maxEProj/minEProj,0.2);
268 
269  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
270  while (E1 <maxEProj*0.9999999){
271  //G4cout<<E1<<'\t'<<E2<<G4endl;
272 
273  int_cross_section +=integral.Simpson(this,
274  &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 5);
275  log_ESec_vector->push_back(std::log(std::min(E2,maxEProj)));
276  log_Prob_vector->push_back(std::log(int_cross_section));
277  E1=E2;
278  E2*=fE;
279 
280  }
281  std::vector< std::vector<G4double>* > res_mat;
282  res_mat.clear();
283  if (int_cross_section >0.) {
284  res_mat.push_back(log_ESec_vector);
285  res_mat.push_back(log_Prob_vector);
286  }
287 
288  return res_mat;
289 }
290 
292 //
294  G4double kinEnergyScatProj,
295  G4double Z,
296  G4double A ,
297  G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
301  kinEnergyScatProjForIntegration = kinEnergyScatProj;
302 
303  //compute the vector of integrated cross sections
304  //-------------------
305 
306  G4double minEProj= GetSecondAdjEnergyMinForScatProjToProjCase(kinEnergyScatProj);
307  G4double maxEProj= GetSecondAdjEnergyMaxForScatProjToProjCase(kinEnergyScatProj);
308  G4double dEmax=maxEProj-kinEnergyScatProj;
309  G4double dEmin=GetLowEnergyLimit();
310  G4double dE1=dEmin;
311  G4double dE2=dEmin;
312 
313 
314  std::vector< double>* log_ESec_vector = new std::vector< double>();
315  std::vector< double>* log_Prob_vector = new std::vector< double>();
316  log_ESec_vector->push_back(std::log(dEmin));
317  log_Prob_vector->push_back(-50.);
318  G4int nbins=std::max( int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
319  G4double fE=std::pow(dEmax/dEmin,1./nbins);
320 
321 
322 
323 
324 
325  G4double int_cross_section=0.;
326 
327  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
328  while (dE1 <dEmax*0.9999999999999){
329  dE2=dE1*fE;
330  int_cross_section +=integral.Simpson(this,
331  &G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 5);
332  //G4cout<<"int_cross_section "<<minEProj+dE1<<'\t'<<int_cross_section<<G4endl;
333  log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj-minEProj)));
334  log_Prob_vector->push_back(std::log(int_cross_section));
335  dE1=dE2;
336 
337  }
338 
339 
340  std::vector< std::vector<G4double> *> res_mat;
341  res_mat.clear();
342  if (int_cross_section >0.) {
343  res_mat.push_back(log_ESec_vector);
344  res_mat.push_back(log_Prob_vector);
345  }
346 
347  return res_mat;
348 }
350 //
352  G4Material* aMaterial,
353  G4double kinEnergyProd,
354  G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
356  SelectedMaterial= aMaterial;
357  kinEnergyProdForIntegration = kinEnergyProd;
358  //compute the vector of integrated cross sections
359  //-------------------
360 
361  G4double minEProj= GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
362  G4double maxEProj= GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
363  G4double E1=minEProj;
364  std::vector< double>* log_ESec_vector = new std::vector< double>();
365  std::vector< double>* log_Prob_vector = new std::vector< double>();
366  log_ESec_vector->clear();
367  log_Prob_vector->clear();
368  log_ESec_vector->push_back(std::log(E1));
369  log_Prob_vector->push_back(-50.);
370 
371  G4double E2=std::pow(10.,double( int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
372  G4double fE=std::pow(10.,1./nbin_pro_decade);
373  G4double int_cross_section=0.;
374 
375  if (std::pow(fE,5.)>(maxEProj/minEProj)) fE = std::pow(maxEProj/minEProj,0.2);
376 
377  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
378  while (E1 <maxEProj*0.9999999){
379 
380  int_cross_section +=integral.Simpson(this,
381  &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 5);
382  log_ESec_vector->push_back(std::log(std::min(E2,maxEProj)));
383  log_Prob_vector->push_back(std::log(int_cross_section));
384  E1=E2;
385  E2*=fE;
386 
387  }
388  std::vector< std::vector<G4double>* > res_mat;
389  res_mat.clear();
390 
391  if (int_cross_section >0.) {
392  res_mat.push_back(log_ESec_vector);
393  res_mat.push_back(log_Prob_vector);
394  }
395 
396 
397 
398  return res_mat;
399 }
400 
402 //
404  G4Material* aMaterial,
405  G4double kinEnergyScatProj,
406  G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
408  SelectedMaterial= aMaterial;
409  kinEnergyScatProjForIntegration = kinEnergyScatProj;
410 
411  //compute the vector of integrated cross sections
412  //-------------------
413 
414  G4double minEProj= GetSecondAdjEnergyMinForScatProjToProjCase(kinEnergyScatProj);
415  G4double maxEProj= GetSecondAdjEnergyMaxForScatProjToProjCase(kinEnergyScatProj);
416 
417 
418  G4double dEmax=maxEProj-kinEnergyScatProj;
419  G4double dEmin=GetLowEnergyLimit();
420  G4double dE1=dEmin;
421  G4double dE2=dEmin;
422 
423 
424  std::vector< double>* log_ESec_vector = new std::vector< double>();
425  std::vector< double>* log_Prob_vector = new std::vector< double>();
426  log_ESec_vector->push_back(std::log(dEmin));
427  log_Prob_vector->push_back(-50.);
428  G4int nbins=std::max( int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
429  G4double fE=std::pow(dEmax/dEmin,1./nbins);
430 
431  G4double int_cross_section=0.;
432 
433  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
434  while (dE1 <dEmax*0.9999999999999){
435  dE2=dE1*fE;
436  int_cross_section +=integral.Simpson(this,
437  &G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 5);
438  log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj-minEProj)));
439  log_Prob_vector->push_back(std::log(int_cross_section));
440  dE1=dE2;
441 
442  }
443 
444 
445 
446 
447 
448  std::vector< std::vector<G4double> *> res_mat;
449  res_mat.clear();
450  if (int_cross_section >0.) {
451  res_mat.push_back(log_ESec_vector);
452  res_mat.push_back(log_Prob_vector);
453  }
454 
455  return res_mat;
456 }
458 //
459 G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex,G4double aPrimEnergy,G4bool IsScatProjToProjCase)
460 {
461 
462 
463  G4AdjointCSMatrix* theMatrix= (*pOnCSMatrixForProdToProjBackwardScattering)[MatrixIndex];
464  if (IsScatProjToProjCase) theMatrix= (*pOnCSMatrixForScatProjToProjBackwardScattering)[MatrixIndex];
465  std::vector< double>* theLogPrimEnergyVector = theMatrix->GetLogPrimEnergyVector();
466 
467  if (theLogPrimEnergyVector->size() ==0){
468  G4cout<<"No data are contained in the given AdjointCSMatrix!"<<G4endl;
469  G4cout<<"The sampling procedure will be stopped."<<G4endl;
470  return 0.;
471 
472  }
473 
475  G4double aLogPrimEnergy = std::log(aPrimEnergy);
476  size_t ind =theInterpolator->FindPositionForLogVector(aLogPrimEnergy,*theLogPrimEnergyVector);
477 
478 
479  G4double aLogPrimEnergy1,aLogPrimEnergy2;
480  G4double aLogCS1,aLogCS2;
481  G4double log01,log02;
482  std::vector< double>* aLogSecondEnergyVector1 =0;
483  std::vector< double>* aLogSecondEnergyVector2 =0;
484  std::vector< double>* aLogProbVector1=0;
485  std::vector< double>* aLogProbVector2=0;
486  std::vector< size_t>* aLogProbVectorIndex1=0;
487  std::vector< size_t>* aLogProbVectorIndex2=0;
488 
489  theMatrix->GetData(ind, aLogPrimEnergy1,aLogCS1,log01, aLogSecondEnergyVector1,aLogProbVector1,aLogProbVectorIndex1);
490  theMatrix->GetData(ind+1, aLogPrimEnergy2,aLogCS2,log02, aLogSecondEnergyVector2,aLogProbVector2,aLogProbVectorIndex2);
491 
492  G4double rand_var = G4UniformRand();
493  G4double log_rand_var= std::log(rand_var);
494  G4double log_Tcut =std::log(currentTcutForDirectSecond);
495  G4double Esec=0;
496  G4double log_dE1,log_dE2;
497  G4double log_rand_var1,log_rand_var2;
498  G4double log_E1,log_E2;
499  log_rand_var1=log_rand_var;
500  log_rand_var2=log_rand_var;
501 
502  G4double Emin=0.;
503  G4double Emax=0.;
504  if (theMatrix->IsScatProjToProjCase()){ //case where Tcut plays a role
507  G4double dE=0;
508  if (Emin < Emax ){
509  if (ApplyCutInRange) {
510  if (second_part_of_same_type && currentTcutForDirectSecond>aPrimEnergy) return aPrimEnergy;
511 
512  log_rand_var1=log_rand_var+theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector1,*aLogProbVector1);
513  log_rand_var2=log_rand_var+theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector2,*aLogProbVector2);
514 
515  }
516  log_dE1 = theInterpolator->Interpolate(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
517  log_dE2 = theInterpolator->Interpolate(log_rand_var2,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
518  dE=std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_dE1,log_dE2));
519  }
520 
521  Esec = aPrimEnergy +dE;
522  Esec=std::max(Esec,Emin);
523  Esec=std::min(Esec,Emax);
524 
525  }
526  else { //Tcut condition is already full-filled
527 
528  log_E1 = theInterpolator->Interpolate(log_rand_var,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
529  log_E2 = theInterpolator->Interpolate(log_rand_var,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
530 
531  Esec = std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_E1,log_E2));
532  Emin=GetSecondAdjEnergyMinForProdToProjCase(aPrimEnergy);
533  Emax=GetSecondAdjEnergyMaxForProdToProjCase(aPrimEnergy);
534  Esec=std::max(Esec,Emin);
535  Esec=std::min(Esec,Emax);
536 
537  }
538 
539  return Esec;
540 
541 
542 
543 
544 
545 }
546 
548 //
550 { SelectCSMatrix(IsScatProjToProjCase);
551  return SampleAdjSecEnergyFromCSMatrix(indexOfUsedCrossSectionMatrix, aPrimEnergy, IsScatProjToProjCase);
552 }
554 //
555 void G4VEmAdjointModel::SelectCSMatrix(G4bool IsScatProjToProjCase)
556 {
559  else if (!UseOnlyOneMatrixForAllElements) { //Select Material
560  std::vector<G4double>* CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
562  if ( !IsScatProjToProjCase) {
563  CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
565  }
566  G4double rand_var= G4UniformRand();
567  G4double SumCS=0.;
568  size_t ind=0;
569  for (size_t i=0;i<CS_Vs_Element->size();i++){
570  SumCS+=(*CS_Vs_Element)[i];
571  if (rand_var<=SumCS/lastCS){
572  ind=i;
573  break;
574  }
575  }
577  }
578 }
580 //
582 {
583  // here we try to use the rejection method
584  //-----------------------------------------
585 
586  const G4int iimax = 1000;
587  G4double E=0;
588  G4double x,xmin,greject,q;
589  if ( IsScatProjToProjCase){
592  xmin=Emin/Emax;
593  G4double grejmax = DiffCrossSectionPerAtomPrimToScatPrim(Emin,prim_energy,1)*prim_energy;
594 
595  G4int ii =0;
596  do {
597  q = G4UniformRand();
598  x = 1./(q*(1./xmin -1.) +1.);
599  E=x*Emax;
600  greject = DiffCrossSectionPerAtomPrimToScatPrim( E,prim_energy ,1)*prim_energy;
601  ++ii;
602  if(ii >= iimax) { break; }
603  }
604  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
605  while( greject < G4UniformRand()*grejmax );
606 
607  }
608  else {
611  xmin=Emin/Emax;
612  G4double grejmax = DiffCrossSectionPerAtomPrimToSecond(Emin,prim_energy,1);
613  G4int ii =0;
614  do {
615  q = G4UniformRand();
616  x = std::pow(xmin, q);
617  E=x*Emax;
618  greject = DiffCrossSectionPerAtomPrimToSecond( E,prim_energy ,1);
619  ++ii;
620  if(ii >= iimax) { break; }
621  }
622  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
623  while( greject < G4UniformRand()*grejmax );
624 
625 
626 
627  }
628 
629  return E;
630 }
631 
633 //
635  G4double old_weight,
636  G4double adjointPrimKinEnergy,
637  G4double projectileKinEnergy,
638  G4bool IsScatProjToProjCase)
639 {
640  G4double new_weight=old_weight;
641  G4double w_corr =1./CS_biasing_factor;
643 
644 
646  if ( !IsScatProjToProjCase) lastCS=lastAdjointCSForProdToProjCase;
647  if ((adjointPrimKinEnergy-preStepEnergy)/preStepEnergy>0.001){ //Is that in all cases needed???
648  G4double post_stepCS=AdjointCrossSection(currentCouple, adjointPrimKinEnergy
649  ,IsScatProjToProjCase );
650  if (post_stepCS>0 && lastCS>0) w_corr*=post_stepCS/lastCS;
651  }
652 
653  new_weight*=w_corr;
654 
655  //G4cout<<"Post step "<<new_weight<<'\t'<<w_corr<<'\t'<<old_weight<<G4endl;
656  new_weight*=projectileKinEnergy/adjointPrimKinEnergy;//This is needed due to the biasing of diff CS
657  //by the factor adjointPrimKinEnergy/projectileKinEnergy
658 
659 
660 
661  fParticleChange->SetParentWeightByProcess(false);
662  fParticleChange->SetSecondaryWeightByProcess(false);
663  fParticleChange->ProposeParentWeight(new_weight);
664 }
666 //
668 { G4double maxEProj= HighEnergyLimit;
669  if (second_part_of_same_type) maxEProj=std::min(kinEnergyScatProj*2.,HighEnergyLimit);
670  return maxEProj;
671 }
673 //
675 { G4double Emin=PrimAdjEnergy;
676  if (ApplyCutInRange) Emin=PrimAdjEnergy+Tcut;
677  return Emin;
678 }
680 //
682 { return HighEnergyLimit;
683 }
685 //
687 { G4double minEProj=PrimAdjEnergy;
688  if (second_part_of_same_type) minEProj=PrimAdjEnergy*2.;
689  return minEProj;
690 }
692 //
694 { if(couple != currentCouple) {
695  currentCouple = const_cast<G4MaterialCutsCouple*> (couple);
696  currentMaterial = const_cast<G4Material*> (couple->GetMaterial());
697  currentCoupleIndex = couple->GetIndex();
699  size_t idx=56;
700  currentTcutForDirectSecond =0.00000000001;
705  if (idx <56){
706  const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
708  }
709  }
710 
711 
712  }
713 }
715 //
717 { HighEnergyLimit=aVal;
719 }
721 //
723 {
724  LowEnergyLimit=aVal;
726 }
728 //
730 {
734  if (theAdjEquivOfDirectPrimPartDef->GetParticleName() =="adj_gamma")
736 }
Float_t x
Definition: compare.C:6
static const G4double dE
G4double ComputeAdjointCS(G4Material *aMaterial, G4VEmAdjointModel *aModel, G4double PrimEnergy, G4double Tcut, G4bool IsScatProjToProjCase, std::vector< G4double > &AdjointCS_for_each_element)
G4double DiffCrossSectionFunction1(G4double kinEnergyProj)
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4VEmAdjointModel(const G4String &nam)
const XML_Char * name
Definition: expat.h:151
G4double lastAdjointCSForScatProjToProjCase
static const G4double Emin
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerAtomForSecond(G4double kinEnergyProd, G4double Z, G4double A=0., G4int nbin_pro_decade=10)
G4double GetLowEnergyLimit()
G4double DiffCrossSectionFunction2(G4double kinEnergyProj)
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy, G4double Tcut=0)
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerAtomForScatProj(G4double kinEnergyProd, G4double Z, G4double A=0., G4int nbin_pro_decade=10)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:707
G4double LinearInterpolation(G4double &x, G4double &x1, G4double &x2, G4double &y1, G4double &y2)
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerVolumeForScatProj(G4Material *aMaterial, G4double kinEnergyProd, G4int nbin_pro_decade=10)
size_t GetIndex() const
Definition: G4Material.hh:262
G4Material * currentMaterial
#define G4endl
Definition: G4ios.hh:61
const G4String & GetParticleName() const
G4double Interpolate(G4double &x, std::vector< G4double > &x_vec, std::vector< G4double > &y_vec, G4String InterPolMethod="Log")
static G4AdjointPositron * AdjointPositron()
static G4AdjointElectron * AdjointElectron()
void SetSecondaryWeightByProcess(G4bool)
G4double InterpolateForLogVector(G4double &x, std::vector< G4double > &x_vec, std::vector< G4double > &y_vec)
G4bool UseOnlyOneMatrixForAllElements
G4double lastAdjointCSForProdToProjCase
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerVolumeForSecond(G4Material *aMaterial, G4double kinEnergyProd, G4int nbin_pro_decade=10)
size_t GetIndex() const
Definition: G4Element.hh:182
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
G4double additional_weight_correction_factor_for_post_step_outside_model
static G4AdjointGamma * AdjointGamma()
G4double currentTcutForDirectSecond
G4MaterialCutsCouple * currentCouple
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4double kinEnergyProjForIntegration
Float_t Z
G4double kinEnergyProdForIntegration
G4ParticleDefinition * theAdjEquivOfDirectSecondPartDef
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
virtual G4double DiffCrossSectionPerAtomPrimToSecond(G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A=0.)
virtual G4double DiffCrossSectionPerVolumePrimToSecond(const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyProd)
G4double GetPostStepWeightCorrection()
void SelectCSMatrix(G4bool IsScatProjToProjCase)
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)
#define G4UniformRand()
Definition: Randomize.hh:53
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
double A(double temperature)
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:357
static G4AdjointInterpolator * GetInstance()
size_t FindPositionForLogVector(G4double &x, std::vector< G4double > &x_vec)
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4bool GetData(unsigned int i, G4double &aPrimEnergy, G4double &aCS, G4double &log0, std::vector< double > *&aLogSecondEnergyVector, std::vector< double > *&aLogProbVector, std::vector< size_t > *&aLogProbVectorIndex)
G4double Simpson(T &typeT, F f, G4double a, G4double b, G4int n)
static G4AdjointCSManager * GetAdjointCSManager()
size_t indexOfUsedCrossSectionMatrix
virtual void CorrectPostStepWeight(G4ParticleChange *fParticleChange, G4double old_weight, G4double adjointPrimKinEnergy, G4double projectileKinEnergy, G4bool IsScatProjToProjCase)
static G4ProductionCutsTable * GetProductionCutsTable()
G4VEmModel * theDirectEMModel
size_t RegisterEmAdjointModel(G4VEmAdjointModel *)
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:252
virtual G4double DiffCrossSectionPerVolumePrimToScatPrim(const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyScatProj)
int G4int
Definition: G4Types.hh:78
void SetParentWeightByProcess(G4bool)
void ProposeParentWeight(G4double finalWeight)
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:714
std::vector< G4double > CS_Vs_ElementForProdToProjCase
G4double DiffCrossSectionPerVolumeFunctionForIntegrationOverEkinProj(G4double EkinProd)
static const G4double Emax
virtual G4double SampleAdjSecEnergyFromDiffCrossSectionPerAtom(G4double prim_energy, G4bool IsScatProjToProjCase)
void SetLowEnergyLimit(G4double aVal)
G4GLOB_DLL std::ostream G4cout
G4double kinEnergyScatProjForIntegration
virtual G4double DiffCrossSectionPerAtomPrimToScatPrim(G4double kinEnergyProj, G4double kinEnergyScatProj, G4double Z, G4double A=0.)
G4double SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex, G4double prim_energy, G4bool IsScatProjToProjCase)
void SetHighEnergyLimit(G4double aVal)
const G4Material * GetMaterial() const
G4ParticleDefinition * theDirectPrimaryPartDef
G4Material * SelectedMaterial
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
const G4Element * GetElement(G4int iel) const
Definition: G4Material.hh:203
std::vector< double > * GetLogPrimEnergyVector()
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
virtual G4double GetAdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
std::vector< G4double > CS_Vs_ElementForScatProjToProjCase
void SetAdjointEquivalentOfDirectPrimaryParticleDefinition(G4ParticleDefinition *aPart)
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
virtual ~G4VEmAdjointModel()