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

Last change on this file since 1036 was 966, checked in by garnier, 17 years ago

fichier ajoutes
File size: 33.5 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#include "G4VEmAdjointModel.hh"
27#include "G4AdjointCSManager.hh"
28
29
30#include "G4Integrator.hh"
31#include "G4TrackStatus.hh"
32#include "G4ParticleChange.hh"
33#include "G4AdjointElectron.hh"
34#include "G4AdjointInterpolator.hh"
35
36////////////////////////////////////////////////////////////////////////////////
37//
38G4VEmAdjointModel::G4VEmAdjointModel(const G4String& nam):
39name(nam)
40// lowLimit(0.1*keV), highLimit(100.0*TeV), fluc(0), name(nam), pParticleChange(0)
41{ G4AdjointCSManager::GetAdjointCSManager()->RegisterEmAdjointModel(this);
42 CorrectWeightMode =true;
43 UseMatrix =true;
44 UseMatrixPerElement = true;
45 ApplyCutInRange = true;
46 ApplyBiasing = true;
47 UseOnlyOneMatrixForAllElements = true;
48 IsIonisation =true;
49 CS_biasing_factor =1.;
50 //ApplyBiasing = false;
51}
52////////////////////////////////////////////////////////////////////////////////
53//
54G4VEmAdjointModel::~G4VEmAdjointModel()
55{;}
56////////////////////////////////////////////////////////////////////////////////
57//
58void G4VEmAdjointModel::SampleSecondaries(const G4Track& aTrack,
59 G4bool IsScatProjToProjCase,
60 G4ParticleChange* fParticleChange)
61{
62
63 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
64 //DefineCurrentMaterial(aTrack->GetMaterialCutsCouple());
65 size_t ind=0;
66 if (!UseMatrixPerElement) ind = currentMaterialIndex;
67 //G4cout<<theAdjointPrimary<<std::endl;
68 else if (!UseOnlyOneMatrixForAllElements) { //Select Material
69 std::vector<double>* CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
70 if ( !IsScatProjToProjCase) CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
71 G4double rand_var= G4UniformRand();
72 G4double SumCS=0.;
73 for (size_t i=0;i<CS_Vs_Element->size();i++){
74 SumCS+=(*CS_Vs_Element)[i];
75 if (rand_var<=SumCS/lastCS){
76 ind=i;
77 break;
78 }
79 }
80 ind = currentMaterial->GetElement(ind)->GetIndex();
81 }
82
83
84
85 //Elastic inverse scattering //not correct in all the cases
86 //---------------------------------------------------------
87 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
88 G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
89 G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
90 //G4cout<<adjointPrimKinEnergy<<std::endl;
91 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
92 return;
93 }
94
95 //Sample secondary energy
96 //-----------------------
97 G4double projectileKinEnergy;
98// if (!IsIonisation ) {
99 projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(ind,
100 adjointPrimKinEnergy,
101 IsScatProjToProjCase);
102 //}
103 /*else {
104 projectileKinEnergy = SampleAdjSecEnergyFromDiffCrossSectionPerAtom(adjointPrimKinEnergy,IsScatProjToProjCase);
105 //G4cout<<projectileKinEnergy<<std::endl;
106 }*/
107 //Weight correction
108 //-----------------------
109
110 CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), adjointPrimKinEnergy,projectileKinEnergy);
111
112
113
114
115
116 //Kinematic
117 //---------
118
119 G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
120 G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
121 G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
122
123
124
125 //Companion
126 //-----------
127 G4double companionM0;
128 companionM0=(adjointPrimTotalEnergy-adjointPrimKinEnergy);
129 if (IsScatProjToProjCase) {
130 companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
131 }
132 G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
133 G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;
134
135
136 //Projectile momentum
137 //--------------------
138 G4double P_parallel = (adjointPrimP*adjointPrimP + projectileP2 - companionP2)/(2.*adjointPrimP);
139 G4double P_perp = std::sqrt( projectileP2 - P_parallel*P_parallel);
140 G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
141 G4double phi =G4UniformRand()*2.*3.1415926;
142 G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
143 projectileMomentum.rotateUz(dir_parallel);
144
145
146
147 if (!IsScatProjToProjCase && CorrectWeightMode){ //kill the primary and add a secondary
148 fParticleChange->ProposeTrackStatus(fStopAndKill);
149 fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
150 //G4cout<<"projectileMomentum "<<projectileMomentum<<std::endl;
151 }
152 else {
153 fParticleChange->ProposeEnergy(projectileKinEnergy);
154 fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
155 }
156}
157////////////////////////////////////////////////////////////////////////////////
158//
159void G4VEmAdjointModel::CorrectPostStepWeight(G4ParticleChange* fParticleChange, G4double old_weight, G4double , G4double )
160{
161 G4double new_weight=old_weight;
162 if (CorrectWeightMode) {
163 G4double w_corr =1./CS_biasing_factor;
164 //G4cout<<w_corr<<std::endl;
165
166 /*G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection(theAdjEquivOfDirectPrimPartDef,
167 theAdjEquivOfDirectSecondPartDef,
168 adjointPrimKinEnergy,projectileKinEnergy,
169 aTrack.GetMaterialCutsCouple());
170 w_corr = projectileKinEnergy;
171 G4double Emin,Emax;
172 if (IsScatProjToProjCase) {
173 Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
174 Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy, currentTcutForDirectSecond);
175
176 }
177 else {
178 Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
179 Emin = GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);
180 }
181 w_corr *=std::log(Emax/Emin)/(Emax-Emin); */
182
183 new_weight*=w_corr;
184 }
185 G4cout<< "new weight"<<new_weight<<std::endl;
186 fParticleChange->SetParentWeightByProcess(false);
187 fParticleChange->SetSecondaryWeightByProcess(false);
188 fParticleChange->ProposeParentWeight(new_weight);
189}
190////////////////////////////////////////////////////////////////////////////////
191//
192G4double G4VEmAdjointModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
193 G4double primEnergy,
194 G4bool IsScatProjToProjCase)
195{
196 DefineCurrentMaterial(aCouple);
197 //G4double fwdCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalForwardCS(G4AdjointElectron::AdjointElectron(),primEnergy,aCouple);
198 //G4double adjCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalAdjointCS(G4AdjointElectron::AdjointElectron(), primEnergy,aCouple);
199 if (IsScatProjToProjCase){
200 lastCS = G4AdjointCSManager::GetAdjointCSManager()->ComputeAdjointCS(currentMaterial,
201 this,
202 primEnergy,
203 currentTcutForDirectSecond,
204 true,
205 CS_Vs_ElementForScatProjToProjCase);
206 /*G4double fwdCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalForwardCS(theAdjEquivOfDirectPrimPartDef,primEnergy,aCouple);
207 G4double adjCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalAdjointCS(theAdjEquivOfDirectPrimPartDef, primEnergy,aCouple);
208 */
209 //if (adjCS >0 )lastCS *=fwdCS/adjCS;
210
211 }
212 else {
213 lastCS = G4AdjointCSManager::GetAdjointCSManager()->ComputeAdjointCS(currentMaterial,
214 this,
215 primEnergy,
216 currentTcutForDirectSecond,
217 false,
218 CS_Vs_ElementForProdToProjCase);
219 /*G4double fwdCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalForwardCS(theAdjEquivOfDirectSecondPartDef,primEnergy,aCouple);
220 G4double adjCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalAdjointCS(theAdjEquivOfDirectSecondPartDef, primEnergy,aCouple);
221 */
222 //if (adjCS >0 )lastCS *=fwdCS/adjCS;
223 //lastCS=0.;
224 }
225
226
227 return lastCS;
228
229}
230////////////////////////////////////////////////////////////////////////////////
231//
232//The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
233G4double G4VEmAdjointModel::DiffCrossSectionPerAtomPrimToSecond(
234 G4double kinEnergyProj,
235 G4double kinEnergyProd,
236 G4double Z,
237 G4double A)
238{
239 G4double dSigmadEprod=0;
240 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
241 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
242
243
244 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile
245 G4double Tmax=kinEnergyProj;
246 if (second_part_of_same_type) Tmax = kinEnergyProj/2.;
247 return Z*DiffCrossSectionMoller(kinEnergyProj,kinEnergyProd);
248 //it could be thta Tmax here should be DBLMAX
249 //Tmax=DBLMAX;
250
251 G4double E1=kinEnergyProd;
252 G4double E2=kinEnergyProd*1.000001;
253 G4double dE=(E2-E1);
254 G4double sigma1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
255 G4double sigma2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
256
257 dSigmadEprod=(sigma1-sigma2)/dE;
258 if (dSigmadEprod>1.) {
259 G4cout<<"sigma1 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma1<<std::endl;
260 G4cout<<"sigma2 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma2<<std::endl;
261 G4cout<<"dsigma "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<dSigmadEprod<<std::endl;
262
263 }
264
265
266
267 }
268
269
270
271
272 return dSigmadEprod;
273
274
275
276}
277//The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
278////////////////////////////////////////////////////////////////////////////////
279//
280G4double G4VEmAdjointModel::DiffCrossSectionPerAtomPrimToScatPrim(
281 G4double kinEnergyProj,
282 G4double kinEnergyScatProj,
283 G4double Z,
284 G4double A)
285{ G4double kinEnergyProd = kinEnergyProj - kinEnergyScatProj;
286 G4double dSigmadEprod;
287 if (kinEnergyProd <=0) dSigmadEprod=0;
288 else dSigmadEprod=DiffCrossSectionPerAtomPrimToSecond(kinEnergyProj,kinEnergyProd,Z,A);
289 return dSigmadEprod;
290
291}
292
293////////////////////////////////////////////////////////////////////////////////
294//
295//The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
296G4double G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToSecond(
297 const G4Material* aMaterial,
298 G4double kinEnergyProj,
299 G4double kinEnergyProd)
300{
301 G4double dSigmadEprod=0;
302 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
303 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
304
305
306 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
307 G4double Tmax=kinEnergyProj;
308 if (second_part_of_same_type) Tmax = kinEnergyProj/2.;
309 //it could be thta Tmax here should be DBLMAX
310 //Tmax=DBLMAX;
311
312 G4double E1=kinEnergyProd;
313
314 G4double E2=kinEnergyProd*1.0001;
315 G4double dE=(E2-E1);
316 G4double sigma1=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E1,E2);
317
318 //G4double sigma2=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E2,1.e50);
319 dSigmadEprod=sigma1/dE;
320 if (dSigmadEprod <0) { //could happen with bremstrahlung dur to suppression effect
321 G4cout<<"Halllllllllllllllllllllllllllllllllllllllllllllllo "<<kinEnergyProj<<'\t'<<E1<<'\t'<<dSigmadEprod<<std::endl;
322 E1=kinEnergyProd;
323 E2=E1*1.1;
324 dE=E2-E1;
325 sigma1=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E1,1.e50);
326 G4double sigma2=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E2,1.e50);
327 dSigmadEprod=(sigma1-sigma2)/dE;
328 G4cout<<dSigmadEprod<<std::endl;
329 }
330
331
332 }
333
334
335
336 return dSigmadEprod;
337
338
339
340}
341//The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
342////////////////////////////////////////////////////////////////////////////////
343//
344G4double G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToScatPrim(
345 const G4Material* aMaterial,
346 G4double kinEnergyProj,
347 G4double kinEnergyScatProj)
348{ G4double kinEnergyProd = kinEnergyProj - kinEnergyScatProj;
349 G4double dSigmadEprod;
350 if (kinEnergyProd <=0) dSigmadEprod=0;
351 else dSigmadEprod=DiffCrossSectionPerVolumePrimToSecond(aMaterial,kinEnergyProj,kinEnergyProd);
352 return dSigmadEprod;
353
354}
355///////////////////////////////////////////////////////////////////////////////////////////////////////////
356//
357G4double G4VEmAdjointModel::DiffCrossSectionFunction1(G4double kinEnergyProj){
358 //return kinEnergyProj*kinEnergyProj;
359 //ApplyBiasing=false;
360 G4double bias_factor = CS_biasing_factor*kinEnergyProdForIntegration/kinEnergyProj;
361 if (!ApplyBiasing) bias_factor =CS_biasing_factor;
362 //G4cout<<bias_factor<<std::endl;
363 if (UseMatrixPerElement ) {
364 return DiffCrossSectionPerAtomPrimToSecond(kinEnergyProj,kinEnergyProdForIntegration,ZSelectedNucleus,ASelectedNucleus)*bias_factor;
365 }
366 else {
367 return DiffCrossSectionPerVolumePrimToSecond(SelectedMaterial,kinEnergyProj,kinEnergyProdForIntegration)*bias_factor;
368 }
369}
370//////////////////////////////////////////////////////////////////////////////
371//
372G4double G4VEmAdjointModel::DiffCrossSectionMoller(G4double kinEnergyProj,G4double kinEnergyProd){
373 G4double electron_mass_c2=0.51099906*MeV;
374 G4double energy = kinEnergyProj + electron_mass_c2;
375 G4double x = kinEnergyProd/kinEnergyProj;
376 G4double gam = energy/electron_mass_c2;
377 G4double gamma2 = gam*gam;
378 G4double beta2 = 1.0 - 1.0/gamma2;
379
380 G4double g = (2.0*gam - 1.0)/gamma2;
381 G4double y = 1.0 - x;
382 G4double fac=twopi_mc2_rcl2/electron_mass_c2;
383 G4double dCS = fac*( 1.-g + ((1.0 - g*x)/(x*x)) + ((1.0 - g*y)/(y*y)))/(beta2*(gam-1));
384 return dCS/kinEnergyProj;
385
386
387
388}
389////////////////////////////////////////////////////////////////////////////////
390//
391G4double G4VEmAdjointModel::DiffCrossSectionFunction2(G4double kinEnergyProj){
392 //return kinEnergyProj*kinEnergyProj;
393 G4double bias_factor = CS_biasing_factor*kinEnergyScatProjForIntegration/kinEnergyProj;
394 //ApplyBiasing=false;
395 if (!ApplyBiasing) bias_factor = CS_biasing_factor;
396 //G4cout<<bias_factor<<std::endl;
397 if (UseMatrixPerElement ) {
398 return DiffCrossSectionPerAtomPrimToScatPrim(kinEnergyProj,kinEnergyScatProjForIntegration,ZSelectedNucleus,ASelectedNucleus)*bias_factor;
399 }
400 else {
401 return DiffCrossSectionPerVolumePrimToScatPrim(SelectedMaterial,kinEnergyProj,kinEnergyScatProjForIntegration)*bias_factor;
402
403 }
404
405}
406////////////////////////////////////////////////////////////////////////////////
407//
408std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForSecond(
409 G4double kinEnergyProd,
410 G4double Z,
411 G4double A ,
412 G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
413{ G4Integrator<G4VEmAdjointModel, G4double(G4VEmAdjointModel::*)(G4double)> integral;
414 ASelectedNucleus= G4int(A);
415 ZSelectedNucleus=G4int(Z);
416 kinEnergyProdForIntegration = kinEnergyProd;
417
418 //compute the vector of integrated cross sections
419 //-------------------
420
421 G4double minEProj= GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
422 G4double maxEProj= GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
423 G4double E1=minEProj;
424 std::vector< G4double >* log_ESec_vector = new std::vector< G4double >();
425 std::vector< G4double >* log_Prob_vector = new std::vector< G4double >();
426 log_ESec_vector->clear();
427 log_Prob_vector->clear();
428 log_ESec_vector->push_back(std::log(E1));
429 log_Prob_vector->push_back(-50.);
430
431 G4double E2=std::pow(10.,G4double( G4int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
432 G4double fE=std::pow(10.,1./nbin_pro_decade);
433 G4double int_cross_section=0.;
434
435 if (std::pow(fE,5.)>(maxEProj/minEProj)) fE = std::pow(maxEProj/minEProj,0.2);
436
437 while (E1 <maxEProj*0.9999999){
438 //G4cout<<E1<<'\t'<<E2<<std::endl;
439
440 int_cross_section +=integral.Simpson(this, &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 10);
441 //G4cout<<"int_cross_section 1 "<<'\t'<<int_cross_section<<std::endl;
442 log_ESec_vector->push_back(std::log(std::min(E2,maxEProj)));
443 log_Prob_vector->push_back(std::log(int_cross_section));
444 E1=E2;
445 E2*=fE;
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}
457
458/////////////////////////////////////////////////////////////////////////////////////
459//
460std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForScatProj(
461 G4double kinEnergyScatProj,
462 G4double Z,
463 G4double A ,
464 G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
465{ G4Integrator<G4VEmAdjointModel, G4double(G4VEmAdjointModel::*)(G4double)> integral;
466 ASelectedNucleus=G4int(A);
467 ZSelectedNucleus=G4int(Z);
468 kinEnergyScatProjForIntegration = kinEnergyScatProj;
469
470 //compute the vector of integrated cross sections
471 //-------------------
472
473 G4double minEProj= GetSecondAdjEnergyMinForScatProjToProjCase(kinEnergyScatProj);
474 G4double maxEProj= GetSecondAdjEnergyMaxForScatProjToProjCase(kinEnergyScatProj);
475 G4double dEmax=maxEProj-kinEnergyScatProj;
476 G4double dEmin=GetLowEnergyLimit();
477 G4double dE1=dEmin;
478 G4double dE2=dEmin;
479
480
481 std::vector< G4double >* log_ESec_vector = new std::vector< G4double >();
482 std::vector< G4double >* log_Prob_vector = new std::vector< G4double >();
483 log_ESec_vector->push_back(std::log(dEmin));
484 log_Prob_vector->push_back(-50.);
485 G4int nbins=std::max( G4int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
486 G4double fE=std::pow(dEmax/dEmin,1./nbins);
487
488 G4double int_cross_section=0.;
489
490 while (dE1 <dEmax*0.9999999999999){
491 dE2=dE1*fE;
492 int_cross_section +=integral.Simpson(this,
493 &G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 20);
494 //G4cout<<"int_cross_section "<<minEProj+dE1<<'\t'<<int_cross_section<<std::endl;
495 log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj)));
496 log_Prob_vector->push_back(std::log(int_cross_section));
497 dE1=dE2;
498
499 }
500 /*G4cout<<"total int_cross_section"<<'\t'<<int_cross_section<<std::endl;
501 G4cout<<"energy "<<kinEnergyScatProj<<std::endl;*/
502
503
504
505
506 std::vector< std::vector<G4double> *> res_mat;
507 res_mat.clear();
508 if (int_cross_section >0.) {
509 res_mat.push_back(log_ESec_vector);
510 res_mat.push_back(log_Prob_vector);
511 }
512
513 return res_mat;
514}
515////////////////////////////////////////////////////////////////////////////////
516//
517std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForSecond(
518 G4Material* aMaterial,
519 G4double kinEnergyProd,
520 G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
521{ G4Integrator<G4VEmAdjointModel, G4double(G4VEmAdjointModel::*)(G4double)> integral;
522 SelectedMaterial= aMaterial;
523 kinEnergyProdForIntegration = kinEnergyProd;
524 //G4cout<<aMaterial->GetName()<<std::endl;
525 //G4cout<<kinEnergyProd/MeV<<std::endl;
526 //compute the vector of integrated cross sections
527 //-------------------
528
529 G4double minEProj= GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
530 G4double maxEProj= GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
531 G4double E1=minEProj;
532 std::vector< G4double >* log_ESec_vector = new std::vector< G4double >();
533 std::vector< G4double >* log_Prob_vector = new std::vector< G4double >();
534 log_ESec_vector->clear();
535 log_Prob_vector->clear();
536 log_ESec_vector->push_back(std::log(E1));
537 log_Prob_vector->push_back(-50.);
538
539 G4double E2=std::pow(10.,G4double( G4int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
540 G4double fE=std::pow(10.,1./nbin_pro_decade);
541 G4double int_cross_section=0.;
542
543 if (std::pow(fE,5.)>(maxEProj/minEProj)) fE = std::pow(maxEProj/minEProj,0.2);
544
545 while (E1 <maxEProj*0.9999999){
546 //G4cout<<E1<<'\t'<<E2<<std::endl;
547
548 int_cross_section +=integral.Simpson(this, &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 10);
549 //G4cout<<"int_cross_section 1 "<<E1<<'\t'<<int_cross_section<<std::endl;
550 log_ESec_vector->push_back(std::log(std::min(E2,maxEProj)));
551 log_Prob_vector->push_back(std::log(int_cross_section));
552 E1=E2;
553 E2*=fE;
554
555 }
556 std::vector< std::vector<G4double>* > res_mat;
557 res_mat.clear();
558
559 //if (int_cross_section >0.) {
560 res_mat.push_back(log_ESec_vector);
561 res_mat.push_back(log_Prob_vector);
562 //}
563
564 return res_mat;
565}
566
567/////////////////////////////////////////////////////////////////////////////////////
568//
569std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForScatProj(
570 G4Material* aMaterial,
571 G4double kinEnergyScatProj,
572 G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
573{ G4Integrator<G4VEmAdjointModel, G4double(G4VEmAdjointModel::*)(G4double)> integral;
574 SelectedMaterial= aMaterial;
575 kinEnergyScatProjForIntegration = kinEnergyScatProj;
576 /*G4cout<<name<<std::endl;
577 G4cout<<aMaterial->GetName()<<std::endl;
578 G4cout<<kinEnergyScatProj/MeV<<std::endl;*/
579 //compute the vector of integrated cross sections
580 //-------------------
581
582 G4double minEProj= GetSecondAdjEnergyMinForScatProjToProjCase(kinEnergyScatProj);
583 G4double maxEProj= GetSecondAdjEnergyMaxForScatProjToProjCase(kinEnergyScatProj);
584
585
586 G4double dEmax=maxEProj-kinEnergyScatProj;
587 G4double dEmin=GetLowEnergyLimit();
588 G4double dE1=dEmin;
589 G4double dE2=dEmin;
590
591
592 std::vector< G4double >* log_ESec_vector = new std::vector< G4double >();
593 std::vector< G4double >* log_Prob_vector = new std::vector< G4double >();
594 log_ESec_vector->push_back(std::log(dEmin));
595 log_Prob_vector->push_back(-50.);
596 G4int nbins=std::max( G4int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
597 G4double fE=std::pow(dEmax/dEmin,1./nbins);
598
599 G4double int_cross_section=0.;
600
601 while (dE1 <dEmax*0.9999999999999){
602 dE2=dE1*fE;
603 int_cross_section +=integral.Simpson(this,
604 &G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 20);
605 //G4cout<<"int_cross_section "<<minEProj+dE1<<'\t'<<int_cross_section<<std::endl;
606 log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj)));
607 log_Prob_vector->push_back(std::log(int_cross_section));
608 dE1=dE2;
609
610 }
611
612
613
614
615
616 std::vector< std::vector<G4double> *> res_mat;
617 res_mat.clear();
618 if (int_cross_section >0.) {
619 res_mat.push_back(log_ESec_vector);
620 res_mat.push_back(log_Prob_vector);
621 }
622
623 return res_mat;
624}
625//////////////////////////////////////////////////////////////////////////////
626//
627G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex,G4double aPrimEnergy,G4bool IsScatProjToProjCase)
628{
629
630
631 G4AdjointCSMatrix* theMatrix= (*pOnCSMatrixForProdToProjBackwardScattering)[MatrixIndex];
632 if (IsScatProjToProjCase) theMatrix= (*pOnCSMatrixForScatProjToProjBackwardScattering)[MatrixIndex];
633 std::vector< G4double >* theLogPrimEnergyVector = theMatrix->GetLogPrimEnergyVector();
634 //G4double dLog = theMatrix->GetDlog();
635
636
637
638 if (theLogPrimEnergyVector->size() ==0){
639 G4cout<<"No data are contained in the given AdjointCSMatrix!"<<std::endl;
640 G4cout<<"The sampling procedure will be stopped."<<std::endl;
641 return 0.;
642
643 }
644
645 G4AdjointInterpolator* theInterpolator=G4AdjointInterpolator::GetInstance();
646 G4double aLogPrimEnergy = std::log(aPrimEnergy);
647 size_t ind =theInterpolator->FindPositionForLogVector(aLogPrimEnergy,*theLogPrimEnergyVector);
648
649
650 G4double aLogPrimEnergy1,aLogPrimEnergy2;
651 G4double aLogCS1,aLogCS2;
652 G4double log01,log02;
653 std::vector< G4double>* aLogSecondEnergyVector1 =0;
654 std::vector< G4double>* aLogSecondEnergyVector2 =0;
655 std::vector< G4double>* aLogProbVector1=0;
656 std::vector< G4double>* aLogProbVector2=0;
657 std::vector< size_t>* aLogProbVectorIndex1=0;
658 std::vector< size_t>* aLogProbVectorIndex2=0;
659
660 theMatrix->GetData(ind, aLogPrimEnergy1,aLogCS1,log01, aLogSecondEnergyVector1,aLogProbVector1,aLogProbVectorIndex1);
661 theMatrix->GetData(ind+1, aLogPrimEnergy2,aLogCS2,log02, aLogSecondEnergyVector2,aLogProbVector2,aLogProbVectorIndex2);
662
663 G4double rand_var = G4UniformRand();
664 G4double log_rand_var= std::log(rand_var);
665 G4double log_Tcut =std::log(currentTcutForDirectSecond);
666 G4double Esec=0;
667 G4double log_dE1,log_dE2;
668 G4double log_rand_var1,log_rand_var2;
669 G4double log_E1,log_E2;
670 log_rand_var1=log_rand_var;
671 log_rand_var2=log_rand_var;
672
673 G4double Emin=0.;
674 G4double Emax=0.;
675 if (theMatrix->IsScatProjToProjCase()){ //case where Tcut plays a role
676 //G4cout<<"Here "<<std::endl;
677 if (ApplyCutInRange) {
678 if (second_part_of_same_type && currentTcutForDirectSecond>aPrimEnergy) return aPrimEnergy;
679 /*if (IsIonisation){
680 G4double inv_Tcut= 1./currentTcutForDirectSecond;
681 G4double inv_dE=inv_Tcut-rand_var*(inv_Tcut-1./aPrimEnergy);
682 Esec= aPrimEnergy+1./inv_dE;
683 //return Esec;
684 G4double dE1=currentTcutForDirectSecond;
685 G4double dE2=currentTcutForDirectSecond*1.00001;
686 G4double dCS1=DiffCrossSectionMoller(aPrimEnergy+dE1,dE1);
687 G4double dCS2=DiffCrossSectionMoller(aPrimEnergy+dE2,dE2);
688 G4double alpha1=std::log(dCS1/dCS2)/std::log(dE1/dE2);
689 G4double a1=dCS1/std::pow(dE1,alpha1);
690 dCS1=DiffCrossSectionMoller(aPrimEnergy+dE1,dE1);
691 dCS2=DiffCrossSectionMoller(aPrimEnergy+dE2,dE2);
692
693 return Esec;
694
695
696
697 dE1=aPrimEnergy/1.00001;
698 dE2=aPrimEnergy;
699 dCS1=DiffCrossSectionMoller(aPrimEnergy+dE1,dE1);
700 dCS2=DiffCrossSectionMoller(aPrimEnergy+dE2,dE2);
701 G4double alpha2=std::log(dCS1/dCS2)/std::log(dE1/dE2);
702 G4double a2=dCS1/std::pow(dE1,alpha1);
703 return Esec;
704
705
706
707
708 }*/
709 log_rand_var1=log_rand_var+theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector1,*aLogProbVector1);
710 log_rand_var2=log_rand_var+theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector2,*aLogProbVector2);
711
712 }
713 log_dE1 = theInterpolator->Interpolate(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
714 log_dE2 = theInterpolator->Interpolate(log_rand_var2,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
715
716 /*log_dE1 = theInterpolator->InterpolateWithIndexVector(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,*aLogProbVectorIndex1,log01,dLog);
717 log_dE2 = theInterpolator->InterpolateWithIndexVector(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,*aLogProbVectorIndex1,log02,dLog);
718 */
719
720
721
722
723
724 Esec = aPrimEnergy +
725 std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_dE1,log_dE2));
726
727 Emin=GetSecondAdjEnergyMinForScatProjToProjCase(aPrimEnergy);
728 Emax=GetSecondAdjEnergyMaxForScatProjToProjCase(aPrimEnergy);
729 Esec=std::max(Esec,Emin);
730 Esec=std::min(Esec,Emax);
731
732
733 //G4cout<<"Esec "<<Esec<<std::endl;
734 //if (Esec > 2.*aPrimEnergy && second_part_of_same_type) Esec = 2.*aPrimEnergy;
735 }
736 else { //Tcut condition is already full-filled
737 /*G4cout<<"Start "<<std::endl;
738 G4cout<<std::exp((*aLogProbVector1)[0])<<std::endl;
739 G4cout<<std::exp((*aLogProbVector2)[0])<<std::endl;*/
740 /*G4double inv_E1= .5/aPrimEnergy;
741
742 G4double inv_E=inv_E1-rand_var*(inv_E1-0.00001);
743 Esec= 1./inv_E;
744 return Esec;*/
745 log_E1 = theInterpolator->Interpolate(log_rand_var,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
746 log_E2 = theInterpolator->Interpolate(log_rand_var,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
747 /*log_E1 = theInterpolator->InterpolateWithIndexVector(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,*aLogProbVectorIndex1,log01,dLog);
748 log_E2 = theInterpolator->InterpolateWithIndexVector(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,*aLogProbVectorIndex1,log02,dLog);
749 */
750
751
752 /*G4cout<<std::exp(log_E1)<<std::endl;
753 G4cout<<std::exp(log_E2)<<std::endl;*/
754
755 Esec = std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_E1,log_E2));
756 Emin=GetSecondAdjEnergyMinForProdToProjCase(aPrimEnergy);
757 Emax=GetSecondAdjEnergyMaxForProdToProjCase(aPrimEnergy);
758 Esec=std::max(Esec,Emin);
759 Esec=std::min(Esec,Emax);
760
761 }
762
763 return Esec;
764
765
766
767
768
769}
770//////////////////////////////////////////////////////////////////////////////
771//
772G4double G4VEmAdjointModel::SampleAdjSecEnergyFromDiffCrossSectionPerAtom(G4double prim_energy,G4bool IsScatProjToProjCase)
773{
774 // here we try to use the rejection method
775 //-----------------------------------------
776
777 G4double E=0;
778 G4double x,xmin,greject,q;
779 if ( IsScatProjToProjCase){
780 G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(prim_energy);
781 G4double Emin= prim_energy+currentTcutForDirectSecond;
782 xmin=Emin/Emax;
783 G4double grejmax = DiffCrossSectionPerAtomPrimToScatPrim(Emin,prim_energy,1)*prim_energy;
784
785 do {
786 q = G4UniformRand();
787 x = 1./(q*(1./xmin -1.) +1.);
788 E=x*Emax;
789 greject = DiffCrossSectionPerAtomPrimToScatPrim( E,prim_energy ,1)*prim_energy;
790
791 }
792
793 while( greject < G4UniformRand()*grejmax );
794
795 }
796 else {
797 G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(prim_energy);
798 G4double Emin= GetSecondAdjEnergyMinForProdToProjCase(prim_energy);;
799 xmin=Emin/Emax;
800 G4double grejmax = DiffCrossSectionPerAtomPrimToSecond(Emin,prim_energy,1);
801 do {
802 q = G4UniformRand();
803 x = std::pow(xmin, q);
804 E=x*Emax;
805 greject = DiffCrossSectionPerAtomPrimToSecond( E,prim_energy ,1);
806
807 }
808
809 while( greject < G4UniformRand()*grejmax );
810
811
812
813 }
814
815 return E;
816
817
818
819
820
821}
822//////////////////////////////////////////////////////////////////////////////
823//
824G4double G4VEmAdjointModel::GetSecondAdjEnergyMaxForScatProjToProjCase(G4double kinEnergyScatProj)
825{ G4double maxEProj= HighEnergyLimit;
826 if (second_part_of_same_type) maxEProj=std::min(kinEnergyScatProj*2.,HighEnergyLimit);
827 return maxEProj;
828}
829//////////////////////////////////////////////////////////////////////////////
830//
831G4double G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy,G4double Tcut)
832{ return PrimAdjEnergy+Tcut;
833}
834//////////////////////////////////////////////////////////////////////////////
835//
836G4double G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProjCase(G4double )
837{ return HighEnergyLimit;
838}
839//////////////////////////////////////////////////////////////////////////////
840//
841G4double G4VEmAdjointModel::GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
842{ G4double minEProj=PrimAdjEnergy;
843 if (second_part_of_same_type) minEProj=PrimAdjEnergy*2.;
844 return minEProj;
845}
846////////////////////////////////////////////////////////////////////////////////////////////
847//
848void G4VEmAdjointModel::DefineCurrentMaterial(const G4MaterialCutsCouple* couple)
849{ if(couple != currentCouple) {
850 currentCouple = const_cast<G4MaterialCutsCouple*> (couple);
851 currentMaterial = const_cast<G4Material*> (couple->GetMaterial());
852 currentCoupleIndex = couple->GetIndex();
853 currentMaterialIndex = currentMaterial->GetIndex();
854 size_t idx=56;
855
856 if (theAdjEquivOfDirectPrimPartDef) {
857 if (theAdjEquivOfDirectPrimPartDef->GetParticleName() == "adj_gamma") idx = 0;
858 else if (theAdjEquivOfDirectPrimPartDef->GetParticleName() == "adj_e-") idx = 1;
859 else if (theAdjEquivOfDirectPrimPartDef->GetParticleName() == "adj_e+") idx = 2;
860 const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
861 currentTcutForDirectPrim=(*aVec)[currentCoupleIndex];
862 }
863
864
865 if (theAdjEquivOfDirectPrimPartDef == theAdjEquivOfDirectSecondPartDef) {
866 currentTcutForDirectSecond = currentTcutForDirectPrim;
867 }
868 else {
869 if (theAdjEquivOfDirectSecondPartDef){
870 if (theAdjEquivOfDirectSecondPartDef->GetParticleName() == "adj_gamma") idx = 0;
871 else if (theAdjEquivOfDirectSecondPartDef->GetParticleName() == "adj_e-") idx = 1;
872 else if (theAdjEquivOfDirectSecondPartDef->GetParticleName() == "adj_e+") idx = 2;
873 const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
874 currentTcutForDirectSecond=(*aVec)[currentCoupleIndex];
875 }
876 }
877 }
878}
Note: See TracBrowser for help on using the repository browser.