Context Navigation

source: trunk/source/processes/electromagnetic/adjoint/src/G4AdjointCSManager.cc @ 1183

Last change on this file since 1183 was 966, checked in by garnier, 15 years ago
fichier ajoutes
File size: 32.1 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 "G4AdjointCSManager.hh"
27	#include "G4AdjointCSMatrix.hh"
28	#include "G4AdjointInterpolator.hh"
29	#include "G4AdjointCSMatrix.hh"
30	#include "G4VEmAdjointModel.hh"
31	#include "G4ElementTable.hh"
32	#include "G4Element.hh"
33	#include "G4ParticleDefinition.hh"
34	#include "G4Element.hh"
35	#include "G4VEmProcess.hh"
36	#include "G4VEnergyLossProcess.hh"
37	#include "G4PhysicsTable.hh"
38	#include "G4PhysicsLogVector.hh"
39	#include "G4PhysicsTableHelper.hh"
40	#include "G4Electron.hh"
41	#include "G4Gamma.hh"
42	#include "G4AdjointElectron.hh"
43	#include "G4AdjointGamma.hh"
44	#include "G4ProductionCutsTable.hh"
45	#include "G4ProductionCutsTable.hh"
46
47
48	G4AdjointCSManager* G4AdjointCSManager::theInstance = 0;
49	///////////////////////////////////////////////////////
50	//
51	G4AdjointCSManager* G4AdjointCSManager::GetAdjointCSManager()
52	{ if(theInstance == 0) {
53	static G4AdjointCSManager ins;
54	theInstance = &ins;
55	}
56	return theInstance;
57	}
58
59	///////////////////////////////////////////////////////
60	//
61	G4AdjointCSManager::G4AdjointCSManager()
62	{ CrossSectionMatrixesAreBuilt=false;
63	theTotalForwardSigmaTableVector.clear();
64	theTotalAdjointSigmaTableVector.clear();
65	listOfForwardEmProcess.clear();
66	listOfForwardEnergyLossProcess.clear();
67	theListOfAdjointParticlesInAction.clear();
68	Tmin=0.1*keV;
69	Tmax=100.*TeV;
70	nbins=240;
71
72	RegisterAdjointParticle(G4AdjointElectron::AdjointElectron());
73	RegisterAdjointParticle(G4AdjointGamma::AdjointGamma());
74
75	verbose = 1;
76
77	consider_continuous_weight_correction =true;
78	consider_poststep_weight_correction =false;
79
80	}
81	///////////////////////////////////////////////////////
82	//
83	G4AdjointCSManager::~G4AdjointCSManager()
84	{;
85	}
86	///////////////////////////////////////////////////////
87	//
88	void G4AdjointCSManager::RegisterEmAdjointModel(G4VEmAdjointModel* aModel)
89	{listOfAdjointEMModel.push_back(aModel);
90	}
91	///////////////////////////////////////////////////////
92	//
93	void G4AdjointCSManager::RegisterEmProcess(G4VEmProcess* aProcess, G4ParticleDefinition* aFwdPartDef)
94	{
95	G4ParticleDefinition* anAdjPartDef = GetAdjointParticleEquivalent(aFwdPartDef);
96	if (anAdjPartDef && aProcess){
97	RegisterAdjointParticle(anAdjPartDef);
98	int index=-1;
99
100	for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
101	if (anAdjPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i;
102	}
103	listOfForwardEmProcess[index]->push_back(aProcess);
104	}
105	}
106	///////////////////////////////////////////////////////
107	//
108	void G4AdjointCSManager::RegisterEnergyLossProcess(G4VEnergyLossProcess* aProcess, G4ParticleDefinition* aFwdPartDef)
109	{
110	G4ParticleDefinition* anAdjPartDef = GetAdjointParticleEquivalent(aFwdPartDef);
111	if (anAdjPartDef && aProcess){
112	RegisterAdjointParticle(anAdjPartDef);
113	int index=-1;
114	for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
115	if (anAdjPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i;
116	}
117	listOfForwardEnergyLossProcess[index]->push_back(aProcess);
118	}
119	}
120	///////////////////////////////////////////////////////
121	//
122	void G4AdjointCSManager::RegisterAdjointParticle(G4ParticleDefinition* aPartDef)
123	{ int index=-1;
124	for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
125	if (aPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i;
126	}
127
128	if (index ==-1){
129	listOfForwardEnergyLossProcess.push_back(new std::vector<G4VEnergyLossProcess*>());
130	theTotalForwardSigmaTableVector.push_back(new G4PhysicsTable);
131	theTotalAdjointSigmaTableVector.push_back(new G4PhysicsTable);
132	listOfForwardEmProcess.push_back(new std::vector<G4VEmProcess*>());
133	theListOfAdjointParticlesInAction.push_back(aPartDef);
134	}
135	}
136	///////////////////////////////////////////////////////
137	//
138	void G4AdjointCSManager::BuildCrossSectionMatrices()
139	{
140	if (CrossSectionMatrixesAreBuilt) return;
141	//Tcut, Tmax
142	//The matrices will be computed probably just once
143	//When Tcut will change some PhysicsTable will be recomputed
144	// for each MaterialCutCouple but not all the matrices
145	//The Tcut defines a lower limit in the energy of the Projectile before the scattering
146	//In the Projectile to Scattered Projectile case we have
147	// E_ScatProj<E_Proj-Tcut
148	//Therefore in the adjoint case we have
149	// Eproj> E_ScatProj+Tcut
150	//This implies that when computing the adjoint CS we should integrate over Epro
151	// from E_ScatProj+Tcut to Emax
152	//In the Projectile to Secondary case Tcut plays a role only in the fact that
153	// Esecond should be greater than Tcut to have the possibility to have any adjoint
154	//process
155	//To avoid to recompute the matrices for all changes of MaterialCutCouple
156	//We propose to compute the matrices only once for the minimum possible Tcut and then
157	//to interpolate the probability for a new Tcut (implemented in G4VAdjointEmModel)
158
159
160	theAdjointCSMatricesForScatProjToProj.clear();
161	theAdjointCSMatricesForProdToProj.clear();
162	const G4ElementTable* theElementTable = G4Element::GetElementTable();
163	const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
164	for (size_t i=0; i<listOfAdjointEMModel.size();i++){
165	G4VEmAdjointModel* aModel =listOfAdjointEMModel[i];
166	G4cout<<"Build adjoint cross section matrices for "<<aModel->GetName()<<std::endl;
167	if (aModel->GetUseMatrix()){
168	std::vector<G4AdjointCSMatrix> aListOfMat1 = new std::vector<G4AdjointCSMatrix*>();
169	std::vector<G4AdjointCSMatrix> aListOfMat2 = new std::vector<G4AdjointCSMatrix*>();
170	aListOfMat1->clear();
171	aListOfMat2->clear();
172	if (aModel->GetUseMatrixPerElement()){
173	if (aModel->GetUseOnlyOneMatrixForAllElements()){
174	std::vector<G4AdjointCSMatrix*>
175	two_matrices=BuildCrossSectionsMatricesForAGivenModelAndElement(aModel,1, 1, 10);
176	aListOfMat1->push_back(two_matrices[0]);
177	aListOfMat2->push_back(two_matrices[1]);
178	}
179	else {
180	for (size_t j=0; j<theElementTable->size();j++){
181	G4Element* anElement=(*theElementTable)[j];
182	G4int Z = G4int(anElement->GetZ());
183	G4int A = G4int(anElement->GetA());
184	std::vector<G4AdjointCSMatrix*>
185	two_matrices=BuildCrossSectionsMatricesForAGivenModelAndElement(aModel,Z, A, 10);
186	aListOfMat1->push_back(two_matrices[0]);
187	aListOfMat2->push_back(two_matrices[1]);
188	}
189	}
190	}
191	else { //Per material case
192	for (size_t j=0; j<theMaterialTable->size();j++){
193	G4Material* aMaterial=(*theMaterialTable)[j];
194	std::vector<G4AdjointCSMatrix*>
195	two_matrices=BuildCrossSectionsMatricesForAGivenModelAndMaterial(aModel,aMaterial, 10);
196	aListOfMat1->push_back(two_matrices[0]);
197	aListOfMat2->push_back(two_matrices[1]);
198	}
199
200	}
201	theAdjointCSMatricesForProdToProj.push_back(*aListOfMat1);
202	theAdjointCSMatricesForScatProjToProj.push_back(*aListOfMat2);
203	aModel->SetCSMatrices(aListOfMat1, aListOfMat2);
204	}
205	else { std::vector<G4AdjointCSMatrix*> two_empty_matrices;
206	theAdjointCSMatricesForProdToProj.push_back(two_empty_matrices);
207	theAdjointCSMatricesForScatProjToProj.push_back(two_empty_matrices);
208
209	}
210	}
211	G4cout<<"All adjoint cross section matrices are built "<<std::endl;
212	CrossSectionMatrixesAreBuilt = true;
213	}
214
215
216	///////////////////////////////////////////////////////
217	//
218	void G4AdjointCSManager::BuildTotalSigmaTables()
219	{
220	const G4ProductionCutsTable* theCoupleTable= G4ProductionCutsTable::GetProductionCutsTable();
221	for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
222	G4ParticleDefinition* thePartDef = theListOfAdjointParticlesInAction[i];
223	theTotalForwardSigmaTableVector[i]->clearAndDestroy();
224	theTotalAdjointSigmaTableVector[i]->clearAndDestroy();
225	for (size_t j=0;j<theCoupleTable->GetTableSize();j++){
226	const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
227
228	//make first the total fwd CS table for FwdProcess
229	G4PhysicsVector* aVector = new G4PhysicsLogVector(Tmin, Tmax, nbins);
230	for(size_t l=0; l<aVector->GetVectorLength(); l++) {
231	G4double totCS=0;
232	G4double e=aVector->GetLowEdgeEnergy(l);
233	for (size_t k=0; k<listOfForwardEmProcess[i]->size(); k++){
234	totCS+=(*listOfForwardEmProcess[i])[k]->GetLambda(e, couple);
235	}
236	for (size_t k=0; k<listOfForwardEnergyLossProcess[i]->size(); k++){
237	totCS+=(*listOfForwardEnergyLossProcess[i])[k]->GetLambda(e, couple);
238	}
239	//G4cout<<totCS<<std::endl;
240	aVector->PutValue(l,totCS);
241
242	}
243	theTotalForwardSigmaTableVector[i]->push_back(aVector);
244
245	G4PhysicsVector* aVector1 = new G4PhysicsLogVector(Tmin, Tmax, nbins);
246	for(size_t l=0; l<aVector->GetVectorLength(); l++) {
247	G4double e=aVector->GetLowEdgeEnergy(l);
248	G4double totCS =ComputeTotalAdjointCS(couple,thePartDef,e);
249	//G4cout<<totCS<<std::endl;
250	aVector1->PutValue(l,totCS);
251
252	}
253	theTotalAdjointSigmaTableVector[i]->push_back(aVector1);
254
255	}
256	}
257
258	}
259	///////////////////////////////////////////////////////
260	//
261	G4double G4AdjointCSManager::GetTotalAdjointCS(G4ParticleDefinition* aPartDef, G4double Ekin,
262	const G4MaterialCutsCouple* aCouple)
263	{ DefineCurrentMaterial(aCouple);
264	int index=-1;
265	for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
266	if (aPartDef == theListOfAdjointParticlesInAction[i]) index=i;
267	}
268	if (index == -1) return 0.;
269
270	G4bool b;
271	return (((*theTotalAdjointSigmaTableVector[index])[currentMatIndex])->GetValue(Ekin, b));
272
273
274
275	}
276	///////////////////////////////////////////////////////
277	//
278	G4double G4AdjointCSManager::GetTotalForwardCS(G4ParticleDefinition* aPartDef, G4double Ekin,
279	const G4MaterialCutsCouple* aCouple)
280	{ DefineCurrentMaterial(aCouple);
281	int index=-1;
282	for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
283	if (aPartDef == theListOfAdjointParticlesInAction[i]) index=i;
284	}
285	if (index == -1) return 0.;
286	G4bool b;
287	return (((*theTotalForwardSigmaTableVector[index])[currentMatIndex])->GetValue(Ekin, b));
288
289
290	}
291	///////////////////////////////////////////////////////
292	//
293	G4double G4AdjointCSManager::GetContinuousWeightCorrection(G4ParticleDefinition* aPartDef, G4double PreStepEkin,G4double AfterStepEkin,
294	const G4MaterialCutsCouple* aCouple, G4double step_length)
295	{ //G4double fwdCS = GetTotalForwardCS(aPartDef, AfterStepEkin,aCouple);
296
297	G4double corr_fac = 1.;
298	if (consider_continuous_weight_correction) {
299
300	G4double adjCS = GetTotalAdjointCS(aPartDef, PreStepEkin,aCouple);
301	G4double PrefwdCS;
302	PrefwdCS = GetTotalForwardCS(aPartDef, PreStepEkin,aCouple);
303	G4double fwdCS = GetTotalForwardCS(aPartDef, (AfterStepEkin+PreStepEkin)/2.,aCouple);
304	G4cout<<adjCS<<'\t'<<fwdCS<<std::endl;
305	//if (aPartDef ==G4AdjointGamma::AdjointGamma()) G4cout<<adjCS<<'\t'<<fwdCS<<std::endl;
306	/if (adjCS >0 ) corr_fac = std::exp((PrefwdCS-fwdCS)step_length);
307	else corr_fac = std::exp(-fwdCSstep_length);/
308	corr_fac =std::exp((adjCS-fwdCS)step_length);
309	corr_fac=std::max(corr_fac,1.e-6);
310	corr_fac *=PreStepEkin/AfterStepEkin;
311
312	}
313	G4cout<<"Cont "<<corr_fac<<std::endl;
314	G4cout<<"Ekin0 "<<PreStepEkin<<std::endl;
315	G4cout<<"Ekin1 "<<AfterStepEkin<<std::endl;
316	G4cout<<"step_length "<<step_length<<std::endl;
317	return corr_fac;
318	}
319	///////////////////////////////////////////////////////
320	//
321	G4double G4AdjointCSManager::GetPostStepWeightCorrection(G4ParticleDefinition* , G4ParticleDefinition* ,
322	G4double ,G4double ,
323	const G4MaterialCutsCouple* )
324	{ G4double corr_fac = 1.;
325	if (consider_poststep_weight_correction) {
326	/*G4double fwdCS = GetTotalForwardCS(aSecondPartDef, EkinPrim,aCouple);
327	G4double adjCS = GetTotalAdjointCS(aPrimPartDef, EkinPrim,aCouple);*/
328	//G4double fwd1CS = GetTotalForwardCS(aPrimPartDef, EkinPrim,aCouple);
329	//if (adjCS>0 && fwd1CS>0) adjCS = fwd1CS;
330	//corr_fac =fwdCS*EkinSecond/adjCS/EkinPrim;
331	//corr_fac = adjCS/fwdCS;
332	}
333	return corr_fac;
334	}
335	///////////////////////////////////////////////////////
336	//
337	double G4AdjointCSManager::ComputeAdjointCS(G4Material* aMaterial,
338	G4VEmAdjointModel* aModel,
339	G4double PrimEnergy,
340	G4double Tcut,
341	G4bool IsScatProjToProjCase,
342	std::vector<double>& CS_Vs_Element)
343	{
344
345	G4bool need_to_compute=false;
346	if ( aMaterial!= lastMaterial \|\| PrimEnergy != lastPrimaryEnergy \|\| Tcut != lastTcut){
347	lastMaterial =aMaterial;
348	lastPrimaryEnergy = PrimEnergy;
349	lastTcut=Tcut;
350	listOfIndexOfAdjointEMModelInAction.clear();
351	listOfIsScatProjToProjCase.clear();
352	lastAdjointCSVsModelsAndElements.clear();
353	need_to_compute=true;
354
355	}
356	size_t ind=0;
357	if (!need_to_compute){
358	need_to_compute=true;
359	for (size_t i=0;i<listOfIndexOfAdjointEMModelInAction.size();i++){
360	size_t ind1=listOfIndexOfAdjointEMModelInAction[i];
361	if (aModel == listOfAdjointEMModel[ind1] && IsScatProjToProjCase == listOfIsScatProjToProjCase[i]){
362	need_to_compute=false;
363	CS_Vs_Element = lastAdjointCSVsModelsAndElements[ind];
364	}
365	ind++;
366	}
367	}
368
369	if (need_to_compute){
370	size_t ind_model=0;
371	for (size_t i=0;i<listOfAdjointEMModel.size();i++){
372	if (aModel == listOfAdjointEMModel[i]){
373	ind_model=i;
374	break;
375	}
376	}
377	G4double Tlow=Tcut;
378	if (!listOfAdjointEMModel[ind_model]->GetApplyCutInRange()) Tlow =listOfAdjointEMModel[ind_model]->GetLowEnergyLimit();
379	listOfIndexOfAdjointEMModelInAction.push_back(ind_model);
380	listOfIsScatProjToProjCase.push_back(IsScatProjToProjCase);
381	CS_Vs_Element.clear();
382	if (!aModel->GetUseMatrix()){
383	return aModel->AdjointCrossSection(currentCouple,PrimEnergy,IsScatProjToProjCase);
384
385
386	}
387	else if (aModel->GetUseMatrixPerElement()){
388	size_t n_el = aMaterial->GetNumberOfElements();
389	if (aModel->GetUseOnlyOneMatrixForAllElements()){
390	G4AdjointCSMatrix* theCSMatrix;
391	if (IsScatProjToProjCase){
392	theCSMatrix=theAdjointCSMatricesForScatProjToProj[ind_model][0];
393	}
394	else theCSMatrix=theAdjointCSMatricesForProdToProj[ind_model][0];
395	G4double CS =0.;
396	if (PrimEnergy > Tlow)
397	CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow);
398	G4double factor=0.;
399	for (size_t i=0;i<n_el;i++){
400	size_t ind_el = aMaterial->GetElement(i)->GetIndex();
401	factor+=aMaterial->GetElement(i)->GetZ()*aMaterial->GetVecNbOfAtomsPerVolume()[i];
402	G4AdjointCSMatrix* theCSMatrix;
403	if (IsScatProjToProjCase){
404	theCSMatrix=theAdjointCSMatricesForScatProjToProj[ind_model][ind_el];
405	}
406	else theCSMatrix=theAdjointCSMatricesForProdToProj[ind_model][ind_el];
407	//G4double CS =0.;
408
409	//G4cout<<CS<<std::endl;
410
411	}
412	CS *=factor;
413	CS_Vs_Element.push_back(CS);
414
415	}
416	else {
417	for (size_t i=0;i<n_el;i++){
418	size_t ind_el = aMaterial->GetElement(i)->GetIndex();
419	//G4cout<<aMaterial->GetName()<<std::endl;
420	G4AdjointCSMatrix* theCSMatrix;
421	if (IsScatProjToProjCase){
422	theCSMatrix=theAdjointCSMatricesForScatProjToProj[ind_model][ind_el];
423	}
424	else theCSMatrix=theAdjointCSMatricesForProdToProj[ind_model][ind_el];
425	G4double CS =0.;
426	if (PrimEnergy > Tlow)
427	CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow);
428	//G4cout<<CS<<std::endl;
429	CS_Vs_Element.push_back(CS*(aMaterial->GetVecNbOfAtomsPerVolume()[i]));
430	}
431	}
432
433	}
434	else {
435	size_t ind_mat = aMaterial->GetIndex();
436	G4AdjointCSMatrix* theCSMatrix;
437	if (IsScatProjToProjCase){
438	theCSMatrix=theAdjointCSMatricesForScatProjToProj[ind_model][ind_mat];
439	}
440	else theCSMatrix=theAdjointCSMatricesForProdToProj[ind_model][ind_mat];
441	G4double CS =0.;
442	if (PrimEnergy > Tlow)
443	CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow);
444	CS_Vs_Element.push_back(CS);
445
446
447	}
448	lastAdjointCSVsModelsAndElements.push_back(CS_Vs_Element);
449
450	}
451
452
453	G4double CS=0;
454	for (size_t i=0;i<CS_Vs_Element.size();i++){
455	CS+=CS_Vs_Element[i];
456	}
457
458	return CS;
459
460
461
462
463
464
465
466
467	}
468	///////////////////////////////////////////////////////
469	//
470	G4Element* G4AdjointCSManager::SampleElementFromCSMatrices(G4Material* aMaterial,
471	G4VEmAdjointModel* aModel,
472	G4double PrimEnergy,
473	G4double Tcut,
474	G4bool IsScatProjToProjCase)
475	{ std::vector<double> CS_Vs_Element;
476	G4double CS = ComputeAdjointCS(aMaterial,aModel,PrimEnergy,Tcut,IsScatProjToProjCase,CS_Vs_Element);
477	G4double rand_var= G4UniformRand();
478	G4double SumCS=0.;
479	size_t ind=0;
480	for (size_t i=0;i<CS_Vs_Element.size();i++){
481	SumCS+=CS_Vs_Element[i];
482	if (rand_var<=SumCS/CS){
483	ind=i;
484	break;
485	}
486	}
487
488	return const_cast<G4Element*>(aMaterial->GetElement(ind));
489
490
491
492	}
493	///////////////////////////////////////////////////////
494	//
495	G4double G4AdjointCSManager::ComputeTotalAdjointCS(const G4MaterialCutsCouple* aCouple,
496	G4ParticleDefinition* aPartDef,
497	G4double Ekin)
498	{
499	G4double TotalCS=0.;
500	// G4ParticleDefinition* theDirPartDef = GetForwardParticleEquivalent(aPartDef);
501	DefineCurrentMaterial(aCouple);
502	/* size_t idx=-1;
503	if (theDirPartDef->GetParticleName() == "gamma") idx = 0;
504	else if (theDirPartDef->GetParticleName() == "e-") idx = 1;
505	else if (theDirPartDef->GetParticleName() == "e+") idx = 2;
506
507	//THe tCut computation is wrong this should be on Tcut per model the secondary determioming the Tcut
508	const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
509	//G4cout<<aVec<<std::endl;
510	G4double Tcut =(aVec)[aCouple->GetIndex()];/
511	//G4cout<<"Tcut "<<Tcut<<std::endl;
512	//G4cout<<(*aVec)[0]<<std::endl;
513	// G4double Tcut =converters[idx]->Convert(Rcut,aCouple->GetMaterial());
514
515
516	std::vector<double> CS_Vs_Element;
517	for (size_t i=0; i<listOfAdjointEMModel.size();i++){
518	/G4ParticleDefinition theDirSecondPartDef =
519	GetForwardParticleEquivalent(listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition());
520
521	*/
522
523
524	G4double Tlow=0;
525	if (!listOfAdjointEMModel[i]->GetApplyCutInRange()) Tlow =listOfAdjointEMModel[i]->GetLowEnergyLimit();
526	else {
527	G4ParticleDefinition* theDirSecondPartDef =
528	GetForwardParticleEquivalent(listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition());
529	G4int idx=-1;
530	if (theDirSecondPartDef->GetParticleName() == "gamma") idx = 0;
531	else if (theDirSecondPartDef->GetParticleName() == "e-") idx = 1;
532	else if (theDirSecondPartDef->GetParticleName() == "e+") idx = 2;
533	const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
534	Tlow =(*aVec)[aCouple->GetIndex()];
535
536
537	}
538
539	if ( Ekin<=listOfAdjointEMModel[i]->GetHighEnergyLimit() && Ekin>=listOfAdjointEMModel[i]->GetLowEnergyLimit()){
540	if (aPartDef == listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectPrimaryParticleDefinition()){
541	//G4cout<<"Yes1 before "<<std::endl;
542	TotalCS += ComputeAdjointCS(currentMaterial,
543	listOfAdjointEMModel[i],
544	Ekin, Tlow,true,CS_Vs_Element);
545	//G4cout<<"Yes1 "<<Ekin<<'\t'<<TotalCS<<std::endl;
546	}
547	if (aPartDef == listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition()){
548	TotalCS += ComputeAdjointCS(currentMaterial,
549	listOfAdjointEMModel[i],
550	Ekin, Tlow,false, CS_Vs_Element);
551
552	//G4cout<<"Yes2 "<<TotalCS<<std::endl;
553	}
554
555	}
556	}
557	return TotalCS;
558
559
560	}
561	///////////////////////////////////////////////////////
562	//
563	std::vector<G4AdjointCSMatrix*>
564	G4AdjointCSManager::BuildCrossSectionsMatricesForAGivenModelAndElement(G4VEmAdjointModel* aModel,G4int Z,G4int A,
565	int nbin_pro_decade)
566	{
567	G4AdjointCSMatrix* theCSMatForProdToProjBackwardScattering = new G4AdjointCSMatrix(false);
568	G4AdjointCSMatrix* theCSMatForScatProjToProjBackwardScattering = new G4AdjointCSMatrix(true);
569
570
571	//make the vector of primary energy of the adjoint particle, could try to make this just once ?
572
573	G4double EkinMin =aModel->GetLowEnergyLimit();
574	G4double EkinMaxForScat =aModel->GetHighEnergyLimit()*0.999;
575	G4double EkinMaxForProd =aModel->GetHighEnergyLimit()*0.999;
576	if (aModel->GetSecondPartOfSameType() )EkinMaxForProd =EkinMaxForProd/2.;
577
578
579
580
581
582
583
584	//Product to projectile backward scattering
585	//-----------------------------------------
586	G4double fE=std::pow(10.,1./nbin_pro_decade);
587	G4double E2=std::pow(10.,G4double( G4int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
588	G4double E1=EkinMin;
589	while (E1 <EkinMaxForProd){
590	E1=std::max(EkinMin,E2);
591	E1=std::min(EkinMaxForProd,E1);
592	std::vector< std::vector< G4double >* > aMat= aModel->ComputeAdjointCrossSectionVectorPerAtomForSecond(E1,Z,A,nbin_pro_decade);
593	if (aMat.size()>=2) {
594	std::vector< G4double >* log_ESecVec=aMat[0];
595	std::vector< G4double >* log_CSVec=aMat[1];
596	G4double log_adjointCS=log_CSVec->back();
597	//normalise CSVec such that it becomes a probability vector
598	/for (size_t j=0;j<log_CSVec->size();j++) (log_CSVec)[j]=(*log_CSVec)[j]-log_adjointCS;
599	(log_CSVec)[0]=-90.;/
600
601
602	for (size_t j=0;j<log_CSVec->size();j++) {
603	//G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<std::endl;
604	if (j==0) (*log_CSVec)[j] = 0.;
605	else (log_CSVec)[j]=std::log(1.-std::exp((log_CSVec)[j]-log_adjointCS));
606	//G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<std::endl;
607	}
608	(log_CSVec)[log_CSVec->size()-1]=(log_CSVec)[log_CSVec->size()-2]-1.;
609	theCSMatForProdToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
610	}
611	E1=E2;
612	E2*=fE;
613	}
614
615	//Scattered projectile to projectile backward scattering
616	//-----------------------------------------
617
618	E2=std::pow(10.,G4double( G4int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
619	E1=EkinMin;
620	while (E1 <EkinMaxForScat){
621	E1=std::max(EkinMin,E2);
622	E1=std::min(EkinMaxForScat,E1);
623	std::vector< std::vector< G4double >* > aMat= aModel->ComputeAdjointCrossSectionVectorPerAtomForScatProj(E1,Z,A,nbin_pro_decade);
624	if (aMat.size()>=2) {
625	std::vector< G4double >* log_ESecVec=aMat[0];
626	std::vector< G4double >* log_CSVec=aMat[1];
627	G4double log_adjointCS=log_CSVec->back();
628	//normalise CSVec such that it becomes a probability vector
629	for (size_t j=0;j<log_CSVec->size();j++) {
630	//G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<std::endl;
631	if (j==0) (*log_CSVec)[j] = 0.;
632	else (log_CSVec)[j]=std::log(1.-std::exp((log_CSVec)[j]-log_adjointCS));
633	//G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<std::endl;
634	}
635	(log_CSVec)[log_CSVec->size()-1]=(log_CSVec)[log_CSVec->size()-2]-1.;
636	theCSMatForScatProjToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
637	}
638	E1=E2;
639	E2*=fE;
640	}
641
642
643
644
645
646
647
648	std::vector<G4AdjointCSMatrix*> res;
649	res.clear();
650
651	res.push_back(theCSMatForProdToProjBackwardScattering);
652	res.push_back(theCSMatForScatProjToProjBackwardScattering);
653
654
655	#ifdef TEST_MODE
656	G4String file_name;
657	std::stringstream astream;
658	G4String str_Z;
659	astream<<Z;
660	astream>>str_Z;
661	theCSMatForProdToProjBackwardScattering->Write(aModel->GetName()+G4String("_CSMat_Z")+str_Z+"_ProdToProj.txt");
662	theCSMatForScatProjToProjBackwardScattering->Write(aModel->GetName()+G4String("_CSMat_Z")+str_Z+"_ScatProjToProj.txt");
663
664	/G4AdjointCSMatrix aMat1 = new G4AdjointCSMatrix(false);
665	G4AdjointCSMatrix* aMat2 = new G4AdjointCSMatrix(true);
666
667	aMat1->Read(G4String("test_Z")+str_Z+"_1.txt");
668	aMat2->Read(G4String("test_Z")+str_Z+"_2.txt");
669	aMat1->Write(G4String("test_Z")+str_Z+"_11.txt");
670	aMat2->Write(G4String("test_Z")+str_Z+"_22.txt"); */
671	#endif
672
673	return res;
674
675
676	}
677	///////////////////////////////////////////////////////
678	//
679	std::vector<G4AdjointCSMatrix*>
680	G4AdjointCSManager::BuildCrossSectionsMatricesForAGivenModelAndMaterial(G4VEmAdjointModel* aModel,
681	G4Material* aMaterial,
682	G4int nbin_pro_decade)
683	{
684	G4AdjointCSMatrix* theCSMatForProdToProjBackwardScattering = new G4AdjointCSMatrix(false);
685	G4AdjointCSMatrix* theCSMatForScatProjToProjBackwardScattering = new G4AdjointCSMatrix(true);
686
687
688	//make the vector of primary energy of the adjoint particle, could try to make this just once ?
689
690	G4double EkinMin =aModel->GetLowEnergyLimit();
691	G4double EkinMaxForScat =aModel->GetHighEnergyLimit()*0.999;
692	G4double EkinMaxForProd =aModel->GetHighEnergyLimit()*0.999;
693	if (aModel->GetSecondPartOfSameType() )EkinMaxForProd =EkinMaxForProd/2.;
694
695
696
697
698
699
700
701	//Product to projectile backward scattering
702	//-----------------------------------------
703	G4double fE=std::pow(10.,1./nbin_pro_decade);
704	G4double E2=std::pow(10.,G4double( G4int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
705	G4double E1=EkinMin;
706	while (E1 <EkinMaxForProd){
707	E1=std::max(EkinMin,E2);
708	E1=std::min(EkinMaxForProd,E1);
709	std::vector< std::vector< G4double >* > aMat= aModel->ComputeAdjointCrossSectionVectorPerVolumeForSecond(aMaterial,E1,nbin_pro_decade);
710	if (aMat.size()>=2) {
711	std::vector< G4double >* log_ESecVec=aMat[0];
712	std::vector< G4double >* log_CSVec=aMat[1];
713	G4double log_adjointCS=log_CSVec->back();
714
715	//normalise CSVec such that it becomes a probability vector
716	for (size_t j=0;j<log_CSVec->size();j++) {
717	//G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<std::endl;
718	if (j==0) (*log_CSVec)[j] = 0.;
719	else (log_CSVec)[j]=std::log(1.-std::exp((log_CSVec)[j]-log_adjointCS));
720	//G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<std::endl;
721	}
722	(log_CSVec)[log_CSVec->size()-1]=(log_CSVec)[log_CSVec->size()-2]-1.;
723	theCSMatForProdToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
724	}
725
726
727
728	E1=E2;
729	E2*=fE;
730	}
731
732	//Scattered projectile to projectile backward scattering
733	//-----------------------------------------
734
735	E2=std::pow(10.,G4double( G4int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
736	E1=EkinMin;
737	while (E1 <EkinMaxForScat){
738	E1=std::max(EkinMin,E2);
739	E1=std::min(EkinMaxForScat,E1);
740	std::vector< std::vector< G4double >* > aMat= aModel->ComputeAdjointCrossSectionVectorPerVolumeForScatProj(aMaterial,E1,nbin_pro_decade);
741	if (aMat.size()>=2) {
742	std::vector< G4double >* log_ESecVec=aMat[0];
743	std::vector< G4double >* log_CSVec=aMat[1];
744	G4double log_adjointCS=log_CSVec->back();
745
746	for (size_t j=0;j<log_CSVec->size();j++) {
747	//G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<std::endl;
748	if (j==0) (*log_CSVec)[j] = 0.;
749	else (log_CSVec)[j]=std::log(1.-std::exp((log_CSVec)[j]-log_adjointCS));
750	//G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<std::endl;
751	}
752	(log_CSVec)[log_CSVec->size()-1]=(log_CSVec)[log_CSVec->size()-2]-1.;
753
754	theCSMatForScatProjToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
755	}
756	E1=E2;
757	E2*=fE;
758	}
759
760
761
762
763
764
765
766	std::vector<G4AdjointCSMatrix*> res;
767	res.clear();
768
769	res.push_back(theCSMatForProdToProjBackwardScattering);
770	res.push_back(theCSMatForScatProjToProjBackwardScattering);
771
772	#ifdef TEST_MODE
773	theCSMatForProdToProjBackwardScattering->Write(aModel->GetName()+"_CSMat_"+aMaterial->GetName()+"_ProdToProj.txt");
774	theCSMatForScatProjToProjBackwardScattering->Write(aModel->GetName()+"_CSMat_"+aMaterial->GetName()+"_ScatProjToProj.txt");
775	#endif
776
777
778	return res;
779
780
781	}
782
783	///////////////////////////////////////////////////////
784	//
785	G4ParticleDefinition* G4AdjointCSManager::GetAdjointParticleEquivalent(G4ParticleDefinition* theFwdPartDef)
786	{
787	if (theFwdPartDef->GetParticleName() == "e-") return G4AdjointElectron::AdjointElectron();
788	if (theFwdPartDef->GetParticleName() == "gamma") return G4AdjointGamma::AdjointGamma();
789	return 0;
790	}
791	///////////////////////////////////////////////////////
792	//
793	G4ParticleDefinition* G4AdjointCSManager::GetForwardParticleEquivalent(G4ParticleDefinition* theAdjPartDef)
794	{
795	if (theAdjPartDef->GetParticleName() == "adj_e-") return G4Electron::Electron();
796	if (theAdjPartDef->GetParticleName() == "adj_gamma") return G4Gamma::Gamma();
797	return 0;
798	}
799	///////////////////////////////////////////////////////
800	//
801	void G4AdjointCSManager::DefineCurrentMaterial(const G4MaterialCutsCouple* couple)
802	{
803	if(couple != currentCouple) {
804	currentCouple = const_cast<G4MaterialCutsCouple*> (couple);
805	currentMaterial = const_cast<G4Material*> (couple->GetMaterial());
806	currentMatIndex = couple->GetIndex();
807	//G4cout<<"Index material "<<currentMatIndex<<std::endl;
808	}
809	}
810
811
812
813	///////////////////////////////////////////////////////
814	//
815	double G4AdjointCSManager::ComputeAdjointCS(G4double aPrimEnergy,G4AdjointCSMatrix*
816	anAdjointCSMatrix,G4double Tcut)
817	{
818	std::vector< G4double > *theLogPrimEnergyVector = anAdjointCSMatrix->GetLogPrimEnergyVector();
819	if (theLogPrimEnergyVector->size() ==0){
820	G4cout<<"No data are contained in the given AdjointCSMatrix!"<<std::endl;
821	G4cout<<"The sampling procedure will be stopped."<<std::endl;
822	return 0.;
823
824	}
825	//G4cout<<"A prim/Tcut "<<aPrimEnergy<<'\t'<<Tcut<<std::endl;
826	G4double log_Tcut = std::log(Tcut);
827	G4double log_E =std::log(aPrimEnergy);
828
829	if (aPrimEnergy <= Tcut \|\| log_E > theLogPrimEnergyVector->back()) return 0.;
830
831
832
833	G4AdjointInterpolator* theInterpolator=G4AdjointInterpolator::GetInstance();
834
835	size_t ind =theInterpolator->FindPositionForLogVector(log_E,*theLogPrimEnergyVector);
836	//G4cout<<"Prim energy "<<(*thePrimEnergyVector)[0]<<std::endl;
837	//G4cout<<"Prim energy[ind]"<<(*thePrimEnergyVector)[ind]<<std::endl;
838	//G4cout<<"Prim energy ind"<<ind<<std::endl;
839
840	G4double aLogPrimEnergy1,aLogPrimEnergy2;
841	G4double aLogCS1,aLogCS2;
842	G4double log01,log02;
843	std::vector< G4double>* aLogSecondEnergyVector1 =0;
844	std::vector< G4double>* aLogSecondEnergyVector2 =0;
845	std::vector< G4double>* aLogProbVector1=0;
846	std::vector< G4double>* aLogProbVector2=0;
847	std::vector< size_t>* aLogProbVectorIndex1=0;
848	std::vector< size_t>* aLogProbVectorIndex2=0;
849
850
851	anAdjointCSMatrix->GetData(ind, aLogPrimEnergy1,aLogCS1,log01, aLogSecondEnergyVector1,aLogProbVector1,aLogProbVectorIndex1);
852	anAdjointCSMatrix->GetData(ind+1, aLogPrimEnergy2,aLogCS2,log02, aLogSecondEnergyVector2,aLogProbVector2,aLogProbVectorIndex2);
853	//G4cout<<"aSecondEnergyVector1.size() "<<aSecondEnergyVector1->size()<<std::endl;
854	//G4cout<<aSecondEnergyVector1<<std::endl;
855	//G4cout<<"aSecondEnergyVector2.size() "<<aSecondEnergyVector2->size()<<std::endl;
856	if (anAdjointCSMatrix->IsScatProjToProjCase()){ //case where the Tcut plays a role
857	G4double log_minimum_prob1, log_minimum_prob2;
858
859	//G4cout<<aSecondEnergyVector1->size()<<std::endl;
860	log_minimum_prob1=theInterpolator->InterpolateForLogVector(log_Tcut,aLogSecondEnergyVector1,aLogProbVector1);
861	log_minimum_prob2=theInterpolator->InterpolateForLogVector(log_Tcut,aLogSecondEnergyVector2,aLogProbVector2);
862	//G4cout<<"minimum_prob1 "<< std::exp(log_minimum_prob1)<<std::endl;
863	//G4cout<<"minimum_prob2 "<< std::exp(log_minimum_prob2)<<std::endl;
864	//G4cout<<"Tcut "<<std::endl;
865	aLogCS1+= log_minimum_prob1;
866	aLogCS2+= log_minimum_prob2;
867	}
868
869	G4double log_adjointCS = theInterpolator->LinearInterpolation(log_E,aLogPrimEnergy1,aLogPrimEnergy2,aLogCS1,aLogCS2);
870	return std::exp(log_adjointCS);
871
872
873	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: