Geant4  10.02.p02
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 93358 2015-10-19 13:41:21Z 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;
51 }
53 //
55 {;}
57 //
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;
68  this,
69  primEnergy,
71  IsScatProjToProjCase,
72  *CS_Vs_Element);
73  if (IsScatProjToProjCase) lastAdjointCSForScatProjToProjCase = lastCS;
75 
76 
77 
78  return lastCS;
79 
80 }
82 //
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
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);
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 //
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
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 //
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 //
202 
203 
204  G4double bias_factor = CS_biasing_factor*kinEnergyProdForIntegration/kinEnergyProj;
205 
206 
207  if (UseMatrixPerElement ) {
209  }
210  else {
212  }
213 }
214 
216 //
218 
219  G4double bias_factor = CS_biasing_factor*kinEnergyScatProjForIntegration/kinEnergyProj;
220  if (UseMatrixPerElement ) {
222  }
223  else {
225 
226  }
227 
228 }
230 //
231 
233 {
235 }
237 //
239  G4double kinEnergyProd,
240  G4double Z,
241  G4double A ,
242  G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
243 {
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 //
293  G4double kinEnergyScatProj,
294  G4double Z,
295  G4double A ,
296  G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
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 //
351  G4Material* aMaterial,
352  G4double kinEnergyProd,
353  G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
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 //
403  G4Material* aMaterial,
404  G4double kinEnergyScatProj,
405  G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
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 
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
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 //
549 { SelectCSMatrix(IsScatProjToProjCase);
550  return SampleAdjSecEnergyFromCSMatrix(indexOfUsedCrossSectionMatrix, aPrimEnergy, IsScatProjToProjCase);
551 }
553 //
554 void G4VEmAdjointModel::SelectCSMatrix(G4bool IsScatProjToProjCase)
555 {
558  else if (!UseOnlyOneMatrixForAllElements) { //Select Material
559  std::vector<G4double>* CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
561  if ( !IsScatProjToProjCase) {
562  CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
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  }
576  }
577 }
579 //
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){
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 {
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 //
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;
642 
643 
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 //
667 { G4double maxEProj= HighEnergyLimit;
668  if (second_part_of_same_type) maxEProj=std::min(kinEnergyScatProj*2.,HighEnergyLimit);
669  return maxEProj;
670 }
672 //
674 { G4double Emin=PrimAdjEnergy;
675  if (ApplyCutInRange) Emin=PrimAdjEnergy+Tcut;
676  return Emin;
677 }
679 //
681 { return HighEnergyLimit;
682 }
684 //
686 { G4double minEProj=PrimAdjEnergy;
687  if (second_part_of_same_type) minEProj=PrimAdjEnergy*2.;
688  return minEProj;
689 }
691 //
693 { if(couple != currentCouple) {
694  currentCouple = const_cast<G4MaterialCutsCouple*> (couple);
695  currentMaterial = const_cast<G4Material*> (couple->GetMaterial());
696  currentCoupleIndex = couple->GetIndex();
698  size_t idx=56;
699  currentTcutForDirectSecond =0.00000000001;
704  if (idx <56){
705  const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
707  }
708  }
709 
710 
711  }
712 }
714 //
716 { HighEnergyLimit=aVal;
718 }
720 //
722 {
723  LowEnergyLimit=aVal;
725 }
727 //
729 {
733  if (theAdjEquivOfDirectPrimPartDef->GetParticleName() =="adj_gamma")
735 }
static G4AdjointGamma * AdjointGamma()
G4double kinEnergyScatProjForIntegration
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:258
virtual ~G4VEmAdjointModel()
virtual G4double DiffCrossSectionPerVolumePrimToSecond(const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyProd)
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)
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 lastAdjointCSForScatProjToProjCase
G4double GetPostStepWeightCorrection()
G4double SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex, G4double prim_energy, G4bool IsScatProjToProjCase)
G4ParticleDefinition * theDirectPrimaryPartDef
size_t GetIndex() const
Definition: G4Material.hh:262
G4String name
Definition: TRTMaterials.hh:40
void ProposeParentWeight(G4double finalWeight)
virtual G4double GetAdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
size_t RegisterEmAdjointModel(G4VEmAdjointModel *)
virtual G4double SampleAdjSecEnergyFromDiffCrossSectionPerAtom(G4double prim_energy, G4bool IsScatProjToProjCase)
G4double DiffCrossSectionPerVolumeFunctionForIntegrationOverEkinProj(G4double EkinProd)
virtual G4double DiffCrossSectionPerVolumePrimToScatPrim(const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyScatProj)
static G4AdjointElectron * AdjointElectron()
virtual G4double DiffCrossSectionPerAtomPrimToSecond(G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A=0.)
G4Material * currentMaterial
const G4Element * GetElement(G4int iel) const
Definition: G4Material.hh:202
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
void SetLowEnergyLimit(G4double aVal)
G4Material * SelectedMaterial
int G4int
Definition: G4Types.hh:78
size_t indexOfUsedCrossSectionMatrix
const G4String & GetParticleName() const
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:725
G4double InterpolateForLogVector(G4double &x, std::vector< G4double > &x_vec, std::vector< G4double > &y_vec)
void SetHighEnergyLimit(G4double aVal)
G4double DiffCrossSectionFunction1(G4double kinEnergyProj)
G4double Interpolate(G4double &x, std::vector< G4double > &x_vec, std::vector< G4double > &y_vec, G4String InterPolMethod="Log")
void SetSecondaryWeightByProcess(G4bool)
void SetParentWeightByProcess(G4bool)
static const G4double dE
G4VEmModel * theDirectEMModel
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
std::vector< G4double > CS_Vs_ElementForProdToProjCase
G4VEmAdjointModel(const G4String &nam)
std::vector< G4double > CS_Vs_ElementForScatProjToProjCase
bool G4bool
Definition: G4Types.hh:79
std::vector< double > * GetLogPrimEnergyVector()
size_t FindPositionForLogVector(G4double &x, std::vector< G4double > &x_vec)
virtual G4double DiffCrossSectionPerAtomPrimToScatPrim(G4double kinEnergyProj, G4double kinEnergyScatProj, G4double Z, G4double A=0.)
size_t GetIndex() const
Definition: G4Element.hh:181
G4double GetLowEnergyLimit()
void SetAdjointEquivalentOfDirectPrimaryParticleDefinition(G4ParticleDefinition *aPart)
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4double DiffCrossSectionFunction2(G4double kinEnergyProj)
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:313
G4MaterialCutsCouple * currentCouple
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
G4bool UseOnlyOneMatrixForAllElements
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerAtomForScatProj(G4double kinEnergyProd, G4double Z, G4double A=0., G4int nbin_pro_decade=10)
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
static G4ProductionCutsTable * GetProductionCutsTable()
T max(const T t1, const T t2)
brief Return the largest of the two arguments
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerAtomForSecond(G4double kinEnergyProd, G4double Z, G4double A=0., G4int nbin_pro_decade=10)
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy, G4double Tcut=0)
G4double LinearInterpolation(G4double &x, G4double &x1, G4double &x2, G4double &y1, G4double &y2)
G4double kinEnergyProjForIntegration
void SelectCSMatrix(G4bool IsScatProjToProjCase)
const G4double x[NPOINTSGL]
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerVolumeForSecond(G4Material *aMaterial, G4double kinEnergyProd, G4int nbin_pro_decade=10)
static const G4double Emin
static G4AdjointInterpolator * GetInstance()
static G4Electron * Electron()
Definition: G4Electron.cc:94
static const G4double Emax
#define G4endl
Definition: G4ios.hh:61
G4ParticleDefinition * theAdjEquivOfDirectSecondPartDef
G4double currentTcutForDirectSecond
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerVolumeForScatProj(G4Material *aMaterial, G4double kinEnergyProd, G4int nbin_pro_decade=10)
G4double lastAdjointCSForProdToProjCase
G4double Simpson(T &typeT, F f, G4double a, G4double b, G4int n)
double G4double
Definition: G4Types.hh:76
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:732
static G4AdjointCSManager * GetAdjointCSManager()
virtual void CorrectPostStepWeight(G4ParticleChange *fParticleChange, G4double old_weight, G4double adjointPrimKinEnergy, G4double projectileKinEnergy, G4bool IsScatProjToProjCase)
G4double kinEnergyProdForIntegration
static G4AdjointPositron * AdjointPositron()
const G4Material * GetMaterial() const
G4double ComputeAdjointCS(G4Material *aMaterial, G4VEmAdjointModel *aModel, G4double PrimEnergy, G4double Tcut, G4bool IsScatProjToProjCase, std::vector< G4double > &AdjointCS_for_each_element)