Geant4  10.00.p01
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 66892 2013-01-17 10:57:59Z gunter $
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  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 //
292  G4double kinEnergyScatProj,
293  G4double Z,
294  G4double A ,
295  G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
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 //
349  G4Material* aMaterial,
350  G4double kinEnergyProd,
351  G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
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 //
400  G4Material* aMaterial,
401  G4double kinEnergyScatProj,
402  G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
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 
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
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 //
545 { SelectCSMatrix(IsScatProjToProjCase);
546  return SampleAdjSecEnergyFromCSMatrix(indexOfUsedCrossSectionMatrix, aPrimEnergy, IsScatProjToProjCase);
547 }
549 //
550 void G4VEmAdjointModel::SelectCSMatrix(G4bool IsScatProjToProjCase)
551 {
554  else if (!UseOnlyOneMatrixForAllElements) { //Select Material
555  std::vector<G4double>* CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
557  if ( !IsScatProjToProjCase) {
558  CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
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  }
572  }
573 }
575 //
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){
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 {
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 //
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;
633 
634 
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 //
658 { G4double maxEProj= HighEnergyLimit;
659  if (second_part_of_same_type) maxEProj=std::min(kinEnergyScatProj*2.,HighEnergyLimit);
660  return maxEProj;
661 }
663 //
665 { G4double Emin=PrimAdjEnergy;
666  if (ApplyCutInRange) Emin=PrimAdjEnergy+Tcut;
667  return Emin;
668 }
670 //
672 { return HighEnergyLimit;
673 }
675 //
677 { G4double minEProj=PrimAdjEnergy;
678  if (second_part_of_same_type) minEProj=PrimAdjEnergy*2.;
679  return minEProj;
680 }
682 //
684 { if(couple != currentCouple) {
685  currentCouple = const_cast<G4MaterialCutsCouple*> (couple);
686  currentMaterial = const_cast<G4Material*> (couple->GetMaterial());
687  currentCoupleIndex = couple->GetIndex();
689  size_t idx=56;
690  currentTcutForDirectSecond =0.00000000001;
695  if (idx <56){
696  const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
698  }
699  }
700 
701 
702  }
703 }
705 //
707 { HighEnergyLimit=aVal;
709 }
711 //
713 {
714  LowEnergyLimit=aVal;
716 }
718 //
720 {
724  if (theAdjEquivOfDirectPrimPartDef->GetParticleName() =="adj_gamma")
726 }
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:245
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:260
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:200
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:683
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:87
G4GLOB_DLL std::ostream G4cout
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
static const G4double A[nN]
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:300
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)
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:690
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)