Context Navigation

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

Last change on this file since 1197 was 1196, checked in by garnier, 16 years ago
update CVS release candidate geant4.9.3.01
File size: 27.2 KB

Rev	Line
[966]	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	//
[1196]	26	// $Id: G4VEmAdjointModel.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
	27	// GEANT4 tag $Name: geant4-09-03-cand-01 $
	28	//
[966]	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"
[1196]	36	#include "G4PhysicsTable.hh"
[966]	37
	38	////////////////////////////////////////////////////////////////////////////////
	39	//
	40	G4VEmAdjointModel::G4VEmAdjointModel(const G4String& nam):
	41	name(nam)
	42	// lowLimit(0.1keV), highLimit(100.0TeV), fluc(0), name(nam), pParticleChange(0)
[1196]	43	{
	44	G4AdjointCSManager::GetAdjointCSManager()->RegisterEmAdjointModel(this);
	45	second_part_of_same_type =false;
	46	theDirectEMModel=0;
[966]	47	}
	48	////////////////////////////////////////////////////////////////////////////////
	49	//
	50	G4VEmAdjointModel::~G4VEmAdjointModel()
	51	{;}
	52	////////////////////////////////////////////////////////////////////////////////
	53	//
	54	G4double G4VEmAdjointModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
	55	G4double primEnergy,
	56	G4bool IsScatProjToProjCase)
	57	{
	58	DefineCurrentMaterial(aCouple);
[1196]	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,
[966]	64	this,
	65	primEnergy,
	66	currentTcutForDirectSecond,
[1196]	67	IsScatProjToProjCase,
	68	*CS_Vs_Element);
	69	if (IsScatProjToProjCase) lastAdjointCSForScatProjToProjCase = lastCS;
	70	else lastAdjointCSForProdToProjCase =lastCS;
	71
[966]	72
	73
	74	return lastCS;
	75
	76	}
	77	////////////////////////////////////////////////////////////////////////////////
	78	//
[1196]	79	//General implementation correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
[966]	80	G4double 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
[1196]	95	G4double E1=kinEnergyProd;
[966]	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	//
	111	G4double 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
	127	G4double 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.;
[1196]	140	G4double E1=kinEnergyProd;
	141	G4double E2=kinEnergyProd*1.0001;
[966]	142	G4double dE=(E2-E1);
	143	G4double sigma1=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E1,E2);
[1196]	144	G4double sigma2=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E2,1.e50);
	145	dSigmadEprod=(sigma1-sigma2)/dE;
[966]	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	//
	155	G4double 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	//
	168	G4double G4VEmAdjointModel::DiffCrossSectionFunction1(G4double kinEnergyProj){
[1196]	169
	170
[966]	171	G4double bias_factor = CS_biasing_factor*kinEnergyProdForIntegration/kinEnergyProj;
[1196]	172
	173
[966]	174	if (UseMatrixPerElement ) {
	175	return DiffCrossSectionPerAtomPrimToSecond(kinEnergyProj,kinEnergyProdForIntegration,ZSelectedNucleus,ASelectedNucleus)*bias_factor;
	176	}
[1196]	177	else {
[966]	178	return DiffCrossSectionPerVolumePrimToSecond(SelectedMaterial,kinEnergyProj,kinEnergyProdForIntegration)*bias_factor;
	179	}
	180	}
	181
	182	////////////////////////////////////////////////////////////////////////////////
	183	//
	184	G4double G4VEmAdjointModel::DiffCrossSectionFunction2(G4double kinEnergyProj){
[1196]	185
	186	G4double bias_factor = CS_biasing_factor*kinEnergyScatProjForIntegration/kinEnergyProj;
[966]	187	if (UseMatrixPerElement ) {
	188	return DiffCrossSectionPerAtomPrimToScatPrim(kinEnergyProj,kinEnergyScatProjForIntegration,ZSelectedNucleus,ASelectedNucleus)*bias_factor;
	189	}
[1196]	190	else {
[966]	191	return DiffCrossSectionPerVolumePrimToScatPrim(SelectedMaterial,kinEnergyProj,kinEnergyScatProjForIntegration)*bias_factor;
	192
	193	}
	194
	195	}
	196	////////////////////////////////////////////////////////////////////////////////
	197	//
[1196]	198
	199	G4double G4VEmAdjointModel::DiffCrossSectionPerVolumeFunctionForIntegrationOverEkinProj(G4double kinEnergyProd)
	200	{
	201	return DiffCrossSectionPerVolumePrimToSecond(SelectedMaterial,kinEnergyProjForIntegration,kinEnergyProd);
	202	}
	203	////////////////////////////////////////////////////////////////////////////////
	204	//
[966]	205	std::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
[1196]	210	{
	211	G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
	212	ASelectedNucleus= int(A);
	213	ZSelectedNucleus=int(Z);
[966]	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;
[1196]	222	std::vector< double>* log_ESec_vector = new std::vector< double>();
	223	std::vector< double>* log_Prob_vector = new std::vector< double>();
[966]	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
[1196]	229	G4double E2=std::pow(10.,double( int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
[966]	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){
[1196]	236	//G4cout<<E1<<'\t'<<E2<<G4endl;
[966]	237
[1196]	238	int_cross_section +=integral.Simpson(this,
	239	&G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 5);
[966]	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	//
	258	std::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
[1196]	263	{ G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
	264	ASelectedNucleus=int(A);
	265	ZSelectedNucleus=int(Z);
[966]	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
[1196]	279	std::vector< double>* log_ESec_vector = new std::vector< double>();
	280	std::vector< double>* log_Prob_vector = new std::vector< double>();
[966]	281	log_ESec_vector->push_back(std::log(dEmin));
	282	log_Prob_vector->push_back(-50.);
[1196]	283	G4int nbins=std::max( int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
[966]	284	G4double fE=std::pow(dEmax/dEmin,1./nbins);
	285
[1196]	286
	287
	288
	289
[966]	290	G4double int_cross_section=0.;
	291
	292	while (dE1 <dEmax*0.9999999999999){
	293	dE2=dE1*fE;
	294	int_cross_section +=integral.Simpson(this,
[1196]	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)));
[966]	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	//
	315	std::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
[1196]	319	{ G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
[966]	320	SelectedMaterial= aMaterial;
	321	kinEnergyProdForIntegration = kinEnergyProd;
[1196]	322	//compute the vector of integrated cross sections
[966]	323	//-------------------
	324
	325	G4double minEProj= GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
	326	G4double maxEProj= GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
	327	G4double E1=minEProj;
[1196]	328	std::vector< double>* log_ESec_vector = new std::vector< double>();
	329	std::vector< double>* log_Prob_vector = new std::vector< double>();
[966]	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
[1196]	335	G4double E2=std::pow(10.,double( int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
[966]	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){
[1196]	342
	343	int_cross_section +=integral.Simpson(this,
	344	&G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 5);
[966]	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
[1196]	354	if (int_cross_section >0.) {
[966]	355	res_mat.push_back(log_ESec_vector);
	356	res_mat.push_back(log_Prob_vector);
[1196]	357	}
[966]	358
[1196]	359
	360
[966]	361	return res_mat;
	362	}
	363
	364	/////////////////////////////////////////////////////////////////////////////////////
	365	//
	366	std::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
[1196]	370	{ G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
[966]	371	SelectedMaterial= aMaterial;
	372	kinEnergyScatProjForIntegration = kinEnergyScatProj;
[1196]	373
[966]	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
[1196]	387	std::vector< double>* log_ESec_vector = new std::vector< double>();
	388	std::vector< double>* log_Prob_vector = new std::vector< double>();
[966]	389	log_ESec_vector->push_back(std::log(dEmin));
	390	log_Prob_vector->push_back(-50.);
[1196]	391	G4int nbins=std::max( int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
[966]	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,
[1196]	399	&G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 5);
	400	log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj-minEProj)));
[966]	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	//
	421	G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex,G4double aPrimEnergy,G4bool IsScatProjToProjCase)
	422	{
	423
	424
	425	G4AdjointCSMatrix* theMatrix= (*pOnCSMatrixForProdToProjBackwardScattering)[MatrixIndex];
	426	if (IsScatProjToProjCase) theMatrix= (*pOnCSMatrixForScatProjToProjBackwardScattering)[MatrixIndex];
[1196]	427	std::vector< double>* theLogPrimEnergyVector = theMatrix->GetLogPrimEnergyVector();
[966]	428
	429	if (theLogPrimEnergyVector->size() ==0){
[1196]	430	G4cout<<"No data are contained in the given AdjointCSMatrix!"<<G4endl;
	431	G4cout<<"The sampling procedure will be stopped."<<G4endl;
[966]	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;
[1196]	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;
[966]	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
[1196]	467	Emin=GetSecondAdjEnergyMinForScatProjToProjCase(aPrimEnergy,currentTcutForDirectSecond);
	468	Emax=GetSecondAdjEnergyMaxForScatProjToProjCase(aPrimEnergy);
	469	G4double dE=0;
	470	if (Emin < Emax ){
	471	if (ApplyCutInRange) {
[966]	472	if (second_part_of_same_type && currentTcutForDirectSecond>aPrimEnergy) return aPrimEnergy;
[1196]	473
[966]	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
[1196]	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	}
[966]	482
[1196]	483	Esec = aPrimEnergy +dE;
[966]	484	Esec=std::max(Esec,Emin);
	485	Esec=std::min(Esec,Emax);
	486
	487	}
	488	else { //Tcut condition is already full-filled
[1196]	489
[966]	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	}
[1196]	508
[966]	509	//////////////////////////////////////////////////////////////////////////////
	510	//
[1196]	511	G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(G4double aPrimEnergy,G4bool IsScatProjToProjCase)
	512	{ SelectCSMatrix(IsScatProjToProjCase);
	513	return SampleAdjSecEnergyFromCSMatrix(indexOfUsedCrossSectionMatrix, aPrimEnergy, IsScatProjToProjCase);
	514	}
	515	//////////////////////////////////////////////////////////////////////////////
	516	//
	517	void 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	//
[966]	543	G4double 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();
[1196]	574	x = pow(xmin, q);
[966]	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;
[1196]	587	}
	588
	589	////////////////////////////////////////////////////////////////////////////////
	590	//
	591	void 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();
[966]	600
[1196]	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);
[966]	621	}
	622	//////////////////////////////////////////////////////////////////////////////
	623	//
	624	G4double 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	//
	631	G4double G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy,G4double Tcut)
[1196]	632	{ G4double Emin=PrimAdjEnergy;
	633	if (ApplyCutInRange) Emin=PrimAdjEnergy+Tcut;
	634	return Emin;
[966]	635	}
	636	//////////////////////////////////////////////////////////////////////////////
	637	//
	638	G4double G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProjCase(G4double )
	639	{ return HighEnergyLimit;
	640	}
	641	//////////////////////////////////////////////////////////////////////////////
	642	//
	643	G4double G4VEmAdjointModel::GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
	644	{ G4double minEProj=PrimAdjEnergy;
	645	if (second_part_of_same_type) minEProj=PrimAdjEnergy*2.;
	646	return minEProj;
	647	}
	648	////////////////////////////////////////////////////////////////////////////////////////////
	649	//
	650	void 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;
[1196]	657	currentTcutForDirectPrim =0.00000000001;
[966]	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;
[1196]	662	if (idx <56){
	663	const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
	664	currentTcutForDirectPrim=(*aVec)[currentCoupleIndex];
	665	}
[966]	666	}
	667
[1196]	668	currentTcutForDirectSecond =0.00000000001;
[966]	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];
[1196]	679	if (idx <56){
	680	const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
	681	currentTcutForDirectPrim=(*aVec)[currentCoupleIndex];
	682	}
	683
	684
[966]	685	}
	686	}
	687	}
	688	}
[1196]	689	////////////////////////////////////////////////////////////////////////////////////////////
	690	//
	691	void G4VEmAdjointModel::SetHighEnergyLimit(G4double aVal)
	692	{ HighEnergyLimit=aVal;
	693	if (theDirectEMModel) theDirectEMModel->SetHighEnergyLimit( aVal);
	694	}
	695	////////////////////////////////////////////////////////////////////////////////////////////
	696	//
	697	void G4VEmAdjointModel::SetLowEnergyLimit(G4double aVal)
	698	{
	699	LowEnergyLimit=aVal;
	700	if (theDirectEMModel) theDirectEMModel->SetLowEnergyLimit( aVal);
	701	}
	702	////////////////////////////////////////////////////////////////////////////////////////////
	703	//
	704	void 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.

Download in other formats: