source: trunk/source/processes/electromagnetic/adjoint/src/G4VEmAdjointModel.cc@ 1340

Last change on this file since 1340 was 1337, checked in by garnier, 15 years ago

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