Context Navigation

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

Last change on this file since 1315 was 1228, checked in by garnier, 15 years ago
update geant4.9.3 tag
File size: 35.7 KB

Line
1	//
2	// ********************************************************************
3	// * License and Disclaimer *
4	// * *
5	// * The Geant4 software is copyright of the Copyright Holders of *
6	// * the Geant4 Collaboration. It is provided under the terms and *
7	// * conditions of the Geant4 Software License, included in the file *
8	// * LICENSE and available at http://cern.ch/geant4/license . These *
9	// * include a list of copyright holders. *
10	// * *
11	// * Neither the authors of this software system, nor their employing *
12	// * institutes,nor the agencies providing financial support for this *
13	// * work make any representation or warranty, express or implied, *
14	// * regarding this software system or assume any liability for its *
15	// * use. Please see the license in the file LICENSE and URL above *
16	// * for the full disclaimer and the limitation of liability. *
17	// * *
18	// * This code implementation is the result of the scientific and *
19	// * technical work of the GEANT4 collaboration. *
20	// * By using, copying, modifying or distributing the software (or *
21	// * any work based on the software) you agree to acknowledge its *
22	// * use in resulting scientific publications, and indicate your *
23	// * acceptance of all terms of the Geant4 Software license. *
24	// ********************************************************************
25	//
26	// $Id: G4AdjointCSManager.cc,v 1.5 2009/11/20 10:31:20 ldesorgh Exp $
27	// GEANT4 tag $Name: geant4-09-03 $
28	//
29	#include "G4AdjointCSManager.hh"
30	#include "G4AdjointCSMatrix.hh"
31	#include "G4AdjointInterpolator.hh"
32	#include "G4AdjointCSMatrix.hh"
33	#include "G4VEmAdjointModel.hh"
34	#include "G4ElementTable.hh"
35	#include "G4Element.hh"
36	#include "G4ParticleDefinition.hh"
37	#include "G4Element.hh"
38	#include "G4VEmProcess.hh"
39	#include "G4VEnergyLossProcess.hh"
40	#include "G4PhysicsTable.hh"
41	#include "G4PhysicsLogVector.hh"
42	#include "G4PhysicsTableHelper.hh"
43	#include "G4Electron.hh"
44	#include "G4Gamma.hh"
45	#include "G4Proton.hh"
46	#include "G4AdjointElectron.hh"
47	#include "G4AdjointGamma.hh"
48	#include "G4AdjointProton.hh"
49	#include "G4ProductionCutsTable.hh"
50	#include "G4ProductionCutsTable.hh"
51	#include <fstream>
52	#include <iomanip>
53
54
55	G4AdjointCSManager* G4AdjointCSManager::theInstance = 0;
56	///////////////////////////////////////////////////////
57	//
58	G4AdjointCSManager* G4AdjointCSManager::GetAdjointCSManager()
59	{ if(theInstance == 0) {
60	static G4AdjointCSManager ins;
61	theInstance = &ins;
62	}
63	return theInstance;
64	}
65
66	///////////////////////////////////////////////////////
67	//
68	G4AdjointCSManager::G4AdjointCSManager()
69	{ CrossSectionMatrixesAreBuilt=false;
70	theTotalForwardSigmaTableVector.clear();
71	theTotalAdjointSigmaTableVector.clear();
72	listOfForwardEmProcess.clear();
73	listOfForwardEnergyLossProcess.clear();
74	theListOfAdjointParticlesInAction.clear();
75	EminForFwdSigmaTables.clear();
76	EminForAdjSigmaTables.clear();
77	EkinofFwdSigmaMax.clear();
78	EkinofAdjSigmaMax.clear();
79	Tmin=0.1*keV;
80	Tmax=100.*TeV;
81	nbins=360; //probably this should be decrease, that was choosen to avoid error in the CS value closed to CS jump.(For example at Tcut)
82
83	RegisterAdjointParticle(G4AdjointElectron::AdjointElectron());
84	RegisterAdjointParticle(G4AdjointGamma::AdjointGamma());
85	RegisterAdjointParticle(G4AdjointProton::AdjointProton());
86
87	verbose = 1;
88
89	lastPartDefForCS =0;
90	LastEkinForCS =0;
91	LastCSCorrectionFactor =1.;
92
93	forward_CS_mode = true;
94
95	currentParticleDef = 0;
96
97	theAdjIon = 0;
98	theFwdIon = 0;
99
100
101	}
102	///////////////////////////////////////////////////////
103	//
104	G4AdjointCSManager::~G4AdjointCSManager()
105	{;
106	}
107	///////////////////////////////////////////////////////
108	//
109	void G4AdjointCSManager::RegisterEmAdjointModel(G4VEmAdjointModel* aModel)
110	{listOfAdjointEMModel.push_back(aModel);
111	}
112	///////////////////////////////////////////////////////
113	//
114	void G4AdjointCSManager::RegisterEmProcess(G4VEmProcess* aProcess, G4ParticleDefinition* aFwdPartDef)
115	{
116	G4ParticleDefinition* anAdjPartDef = GetAdjointParticleEquivalent(aFwdPartDef);
117	if (anAdjPartDef && aProcess){
118	RegisterAdjointParticle(anAdjPartDef);
119	G4int index=-1;
120
121	for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
122	if (anAdjPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i;
123	}
124	listOfForwardEmProcess[index]->push_back(aProcess);
125	}
126	}
127	///////////////////////////////////////////////////////
128	//
129	void G4AdjointCSManager::RegisterEnergyLossProcess(G4VEnergyLossProcess* aProcess, G4ParticleDefinition* aFwdPartDef)
130	{
131	G4ParticleDefinition* anAdjPartDef = GetAdjointParticleEquivalent(aFwdPartDef);
132	if (anAdjPartDef && aProcess){
133	RegisterAdjointParticle(anAdjPartDef);
134	G4int index=-1;
135	for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
136	if (anAdjPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i;
137	}
138	listOfForwardEnergyLossProcess[index]->push_back(aProcess);
139	}
140	}
141	///////////////////////////////////////////////////////
142	//
143	void G4AdjointCSManager::RegisterAdjointParticle(G4ParticleDefinition* aPartDef)
144	{ G4int index=-1;
145	for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
146	if (aPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i;
147	}
148
149	if (index ==-1){
150	listOfForwardEnergyLossProcess.push_back(new std::vector<G4VEnergyLossProcess*>());
151	theTotalForwardSigmaTableVector.push_back(new G4PhysicsTable);
152	theTotalAdjointSigmaTableVector.push_back(new G4PhysicsTable);
153	listOfForwardEmProcess.push_back(new std::vector<G4VEmProcess*>());
154	theListOfAdjointParticlesInAction.push_back(aPartDef);
155	EminForFwdSigmaTables.push_back(std::vector<G4double> ());
156	EminForAdjSigmaTables.push_back(std::vector<G4double> ());
157	EkinofFwdSigmaMax.push_back(std::vector<G4double> ());
158	EkinofAdjSigmaMax.push_back(std::vector<G4double> ());
159
160	}
161	}
162	///////////////////////////////////////////////////////
163	//
164	void G4AdjointCSManager::BuildCrossSectionMatrices()
165	{
166	if (CrossSectionMatrixesAreBuilt) return;
167	//Tcut, Tmax
168	//The matrices will be computed probably just once
169	//When Tcut will change some PhysicsTable will be recomputed
170	// for each MaterialCutCouple but not all the matrices
171	//The Tcut defines a lower limit in the energy of the Projectile before the scattering
172	//In the Projectile to Scattered Projectile case we have
173	// E_ScatProj<E_Proj-Tcut
174	//Therefore in the adjoint case we have
175	// Eproj> E_ScatProj+Tcut
176	//This implies that when computing the adjoint CS we should integrate over Epro
177	// from E_ScatProj+Tcut to Emax
178	//In the Projectile to Secondary case Tcut plays a role only in the fact that
179	// Esecond should be greater than Tcut to have the possibility to have any adjoint
180	//process
181	//To avoid to recompute the matrices for all changes of MaterialCutCouple
182	//We propose to compute the matrices only once for the minimum possible Tcut and then
183	//to interpolate the probability for a new Tcut (implemented in G4VAdjointEmModel)
184
185
186	theAdjointCSMatricesForScatProjToProj.clear();
187	theAdjointCSMatricesForProdToProj.clear();
188	const G4ElementTable* theElementTable = G4Element::GetElementTable();
189	const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
190
191	G4cout<<"========== Computation of cross section matrices for adjoint models =========="<<G4endl;
192	for (size_t i=0; i<listOfAdjointEMModel.size();i++){
193	G4VEmAdjointModel* aModel =listOfAdjointEMModel[i];
194	G4cout<<"Build adjoint cross section matrices for "<<aModel->GetName()<<G4endl;
195	if (aModel->GetUseMatrix()){
196	std::vector<G4AdjointCSMatrix> aListOfMat1 = new std::vector<G4AdjointCSMatrix*>();
197	std::vector<G4AdjointCSMatrix> aListOfMat2 = new std::vector<G4AdjointCSMatrix*>();
198	aListOfMat1->clear();
199	aListOfMat2->clear();
200	if (aModel->GetUseMatrixPerElement()){
201	if (aModel->GetUseOnlyOneMatrixForAllElements()){
202	std::vector<G4AdjointCSMatrix*>
203	two_matrices=BuildCrossSectionsMatricesForAGivenModelAndElement(aModel,1, 1, 80);
204	aListOfMat1->push_back(two_matrices[0]);
205	aListOfMat2->push_back(two_matrices[1]);
206	}
207	else {
208	for (size_t j=0; j<theElementTable->size();j++){
209	G4Element* anElement=(*theElementTable)[j];
210	G4int Z = int(anElement->GetZ());
211	G4int A = int(anElement->GetA());
212	std::vector<G4AdjointCSMatrix*>
213	two_matrices=BuildCrossSectionsMatricesForAGivenModelAndElement(aModel,Z, A, 40);
214	aListOfMat1->push_back(two_matrices[0]);
215	aListOfMat2->push_back(two_matrices[1]);
216	}
217	}
218	}
219	else { //Per material case
220	for (size_t j=0; j<theMaterialTable->size();j++){
221	G4Material* aMaterial=(*theMaterialTable)[j];
222	std::vector<G4AdjointCSMatrix*>
223	two_matrices=BuildCrossSectionsMatricesForAGivenModelAndMaterial(aModel,aMaterial, 40);
224	aListOfMat1->push_back(two_matrices[0]);
225	aListOfMat2->push_back(two_matrices[1]);
226	}
227
228	}
229	theAdjointCSMatricesForProdToProj.push_back(*aListOfMat1);
230	theAdjointCSMatricesForScatProjToProj.push_back(*aListOfMat2);
231	aModel->SetCSMatrices(aListOfMat1, aListOfMat2);
232	}
233	else { G4cout<<"The model "<<aModel->GetName()<<" does not use cross section matrices"<<G4endl;
234	std::vector<G4AdjointCSMatrix*> two_empty_matrices;
235	theAdjointCSMatricesForProdToProj.push_back(two_empty_matrices);
236	theAdjointCSMatricesForScatProjToProj.push_back(two_empty_matrices);
237
238	}
239	}
240	G4cout<<" All adjoint cross section matrices are computed!"<<G4endl;
241	G4cout<<"======================================================================"<<G4endl;
242
243	CrossSectionMatrixesAreBuilt = true;
244
245
246	}
247
248
249	///////////////////////////////////////////////////////
250	//
251	void G4AdjointCSManager::BuildTotalSigmaTables()
252	{
253	const G4ProductionCutsTable* theCoupleTable= G4ProductionCutsTable::GetProductionCutsTable();
254	for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
255	G4ParticleDefinition* thePartDef = theListOfAdjointParticlesInAction[i];
256	DefineCurrentParticle(thePartDef);
257	theTotalForwardSigmaTableVector[i]->clearAndDestroy();
258	theTotalAdjointSigmaTableVector[i]->clearAndDestroy();
259	EminForFwdSigmaTables[i].clear();
260	EminForAdjSigmaTables[i].clear();
261	EkinofFwdSigmaMax[i].clear();
262	EkinofAdjSigmaMax[i].clear();
263	//G4cout<<thePartDef->GetParticleName();
264
265	for (size_t j=0;j<theCoupleTable->GetTableSize();j++){
266	const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
267
268	/*
269	G4String file_name1=couple->GetMaterial()->GetName()+"_"+thePartDef->GetParticleName()+"_adj_totCS.txt";
270	G4String file_name2=couple->GetMaterial()->GetName()+"_"+thePartDef->GetParticleName()+"_fwd_totCS.txt";
271
272	std::fstream FileOutputAdjCS(file_name1, std::ios::out);
273	std::fstream FileOutputFwdCS(file_name2, std::ios::out);
274
275
276
277	FileOutputAdjCS<<std::setiosflags(std::ios::scientific);
278	FileOutputAdjCS<<std::setprecision(6);
279	FileOutputFwdCS<<std::setiosflags(std::ios::scientific);
280	FileOutputFwdCS<<std::setprecision(6);
281	*/
282
283
284	//make first the total fwd CS table for FwdProcess
285	G4PhysicsVector* aVector = new G4PhysicsLogVector(Tmin, Tmax, nbins);
286	G4bool Emin_found=false;
287	size_t ind=0;
288	G4double sigma_max =0.;
289	G4double e_sigma_max =0.;
290	for(size_t l=0; l<aVector->GetVectorLength(); l++) {
291	G4double totCS=0.;
292	G4double e=aVector->GetLowEdgeEnergy(l);
293	for (size_t k=0; k<listOfForwardEmProcess[i]->size(); k++){
294	totCS+=(*listOfForwardEmProcess[i])[k]->GetLambda(e, couple);
295	}
296	for (size_t k=0; k<listOfForwardEnergyLossProcess[i]->size(); k++){
297	if (thePartDef == theAdjIon) { // e is considered already as the scaled energy
298	size_t mat_index = couple->GetIndex();
299	G4VEmModel* currentModel = (*listOfForwardEnergyLossProcess[i])[k]->SelectModelForMaterial(e,mat_index);
300	G4double chargeSqRatio = currentModel->GetChargeSquareRatio(theFwdIon,couple->GetMaterial(),e/massRatio);
301	(*listOfForwardEnergyLossProcess[i])[k]->SetDynamicMassCharge(massRatio,chargeSqRatio);
302	}
303	G4double e1=e/massRatio;
304	totCS+=(*listOfForwardEnergyLossProcess[i])[k]->GetLambda(e1, couple);
305	}
306	aVector->PutValue(l,totCS);
307	if (totCS>sigma_max){
308	sigma_max=totCS;
309	e_sigma_max = e;
310
311	}
312	//FileOutputFwdCS<<e<<'\t'<<totCS<<G4endl;
313
314	if (totCS>0 && !Emin_found) {
315	EminForFwdSigmaTables[i].push_back(e);
316	Emin_found=true;
317	}
318
319
320	}
321	//FileOutputFwdCS.close();
322
323	EkinofFwdSigmaMax[i].push_back(e_sigma_max);
324
325
326	if(!Emin_found) EminForFwdSigmaTables[i].push_back(Tmax);
327
328	theTotalForwardSigmaTableVector[i]->push_back(aVector);
329
330
331	Emin_found=false;
332	sigma_max=0;
333	e_sigma_max =0.;
334	ind=0;
335	G4PhysicsVector* aVector1 = new G4PhysicsLogVector(Tmin, Tmax, nbins);
336	for(size_t l=0; l<aVector->GetVectorLength(); l++) {
337	G4double e=aVector->GetLowEdgeEnergy(l);
338	G4double totCS =ComputeTotalAdjointCS(couple,thePartDef,e*0.9999999/massRatio); //massRatio needed for ions
339	aVector1->PutValue(l,totCS);
340	if (totCS>sigma_max){
341	sigma_max=totCS;
342	e_sigma_max = e;
343
344	}
345	//FileOutputAdjCS<<e<<'\t'<<totCS<<G4endl;
346	if (totCS>0 && !Emin_found) {
347	EminForAdjSigmaTables[i].push_back(e);
348	Emin_found=true;
349	}
350
351	}
352	//FileOutputAdjCS.close();
353	EkinofAdjSigmaMax[i].push_back(e_sigma_max);
354	if(!Emin_found) EminForAdjSigmaTables[i].push_back(Tmax);
355
356	theTotalAdjointSigmaTableVector[i]->push_back(aVector1);
357
358	}
359	}
360
361	}
362	///////////////////////////////////////////////////////
363	//
364	G4double G4AdjointCSManager::GetTotalAdjointCS(G4ParticleDefinition* aPartDef, G4double Ekin,
365	const G4MaterialCutsCouple* aCouple)
366	{ DefineCurrentMaterial(aCouple);
367	DefineCurrentParticle(aPartDef);
368	G4bool b;
369	return (((theTotalAdjointSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(EkinmassRatio, b));
370
371
372
373	}
374	///////////////////////////////////////////////////////
375	//
376	G4double G4AdjointCSManager::GetTotalForwardCS(G4ParticleDefinition* aPartDef, G4double Ekin,
377	const G4MaterialCutsCouple* aCouple)
378	{ DefineCurrentMaterial(aCouple);
379	DefineCurrentParticle(aPartDef);
380	G4bool b;
381	return (((theTotalForwardSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(EkinmassRatio, b));
382
383
384	}
385
386	///////////////////////////////////////////////////////
387	//
388	void G4AdjointCSManager::GetEminForTotalCS(G4ParticleDefinition* aPartDef,
389	const G4MaterialCutsCouple* aCouple, G4double& emin_adj, G4double& emin_fwd)
390	{ DefineCurrentMaterial(aCouple);
391	DefineCurrentParticle(aPartDef);
392	emin_adj = EminForAdjSigmaTables[currentParticleIndex][currentMatIndex]/massRatio;
393	emin_fwd = EminForFwdSigmaTables[currentParticleIndex][currentMatIndex]/massRatio;
394
395
396
397	}
398	///////////////////////////////////////////////////////
399	//
400	void G4AdjointCSManager::GetMaxFwdTotalCS(G4ParticleDefinition* aPartDef,
401	const G4MaterialCutsCouple* aCouple, G4double& e_sigma_max, G4double& sigma_max)
402	{ DefineCurrentMaterial(aCouple);
403	DefineCurrentParticle(aPartDef);
404	e_sigma_max = EkinofFwdSigmaMax[currentParticleIndex][currentMatIndex];
405	G4bool b;
406	sigma_max =((*theTotalForwardSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(e_sigma_max, b);
407	e_sigma_max/=massRatio;
408
409
410	}
411	///////////////////////////////////////////////////////
412	//
413	void G4AdjointCSManager::GetMaxAdjTotalCS(G4ParticleDefinition* aPartDef,
414	const G4MaterialCutsCouple* aCouple, G4double& e_sigma_max, G4double& sigma_max)
415	{ DefineCurrentMaterial(aCouple);
416	DefineCurrentParticle(aPartDef);
417	e_sigma_max = EkinofAdjSigmaMax[currentParticleIndex][currentMatIndex];
418	G4bool b;
419	sigma_max =((*theTotalAdjointSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(e_sigma_max, b);
420	e_sigma_max/=massRatio;
421
422
423	}
424	///////////////////////////////////////////////////////
425	//
426	G4double G4AdjointCSManager::GetCrossSectionCorrection(G4ParticleDefinition* aPartDef,G4double PreStepEkin,const G4MaterialCutsCouple* aCouple, G4bool& fwd_is_used,
427	G4double& fwd_TotCS)
428	{ G4double corr_fac = 1.;
429	if (forward_CS_mode) {
430	fwd_TotCS=PrefwdCS;
431	if (LastEkinForCS != PreStepEkin \|\| aPartDef != lastPartDefForCS \|\| aCouple!=currentCouple) {
432	DefineCurrentMaterial(aCouple);
433	PreadjCS = GetTotalAdjointCS(aPartDef, PreStepEkin,aCouple);
434	PrefwdCS = GetTotalForwardCS(aPartDef, PreStepEkin,aCouple);
435	LastEkinForCS = PreStepEkin;
436	lastPartDefForCS = aPartDef;
437	if (PrefwdCS >0. && PreadjCS >0.) {
438	forward_CS_is_used = true;
439	LastCSCorrectionFactor = PrefwdCS/PreadjCS;
440	}
441	else {
442	forward_CS_is_used = false;
443	LastCSCorrectionFactor = 1.;
444
445	}
446
447	}
448	corr_fac =LastCSCorrectionFactor;
449
450
451
452	}
453	else {
454	forward_CS_is_used = false;
455	LastCSCorrectionFactor = 1.;
456	}
457	fwd_TotCS=PrefwdCS;
458	fwd_is_used = forward_CS_is_used;
459	return corr_fac;
460	}
461	///////////////////////////////////////////////////////
462	//
463	G4double G4AdjointCSManager::GetContinuousWeightCorrection(G4ParticleDefinition* aPartDef, G4double PreStepEkin,G4double AfterStepEkin,
464	const G4MaterialCutsCouple* aCouple, G4double step_length)
465	{ G4double corr_fac = 1.;
466	//return corr_fac;
467	//G4double after_adjCS = GetTotalAdjointCS(aPartDef, AfterStepEkin,aCouple);
468	G4double after_fwdCS = GetTotalForwardCS(aPartDef, AfterStepEkin,aCouple);
469	G4double pre_adjCS = GetTotalAdjointCS(aPartDef, PreStepEkin,aCouple);
470	if (!forward_CS_is_used \|\| pre_adjCS ==0. \|\| after_fwdCS==0.) {
471	forward_CS_is_used=false;
472	G4double pre_fwdCS = GetTotalForwardCS(aPartDef, PreStepEkin,aCouple);
473	corr_fac =std::exp((pre_adjCS-pre_fwdCS)step_length);
474	LastCSCorrectionFactor = 1.;
475	}
476	else {
477	LastCSCorrectionFactor = after_fwdCS/pre_adjCS;
478	}
479
480
481
482	return corr_fac;
483	}
484	///////////////////////////////////////////////////////
485	//
486	G4double G4AdjointCSManager::GetPostStepWeightCorrection( )
487	{//return 1.;
488	return 1./LastCSCorrectionFactor;
489
490	}
491	///////////////////////////////////////////////////////
492	//
493	G4double G4AdjointCSManager::ComputeAdjointCS(G4Material* aMaterial,
494	G4VEmAdjointModel* aModel,
495	G4double PrimEnergy,
496	G4double Tcut,
497	G4bool IsScatProjToProjCase,
498	std::vector<G4double>& CS_Vs_Element)
499	{
500
501	G4double EminSec=0;
502	G4double EmaxSec=0;
503
504	if (IsScatProjToProjCase){
505	EminSec= aModel->GetSecondAdjEnergyMinForScatProjToProjCase(PrimEnergy,Tcut);
506	EmaxSec= aModel->GetSecondAdjEnergyMaxForScatProjToProjCase(PrimEnergy);
507	}
508	else if (PrimEnergy > Tcut \|\| !aModel->GetApplyCutInRange()) {
509	EminSec= aModel->GetSecondAdjEnergyMinForProdToProjCase(PrimEnergy);
510	EmaxSec= aModel->GetSecondAdjEnergyMaxForProdToProjCase(PrimEnergy);
511	}
512	if (EminSec >= EmaxSec) return 0.;
513
514
515	G4bool need_to_compute=false;
516	if ( aMaterial!= lastMaterial \|\| PrimEnergy != lastPrimaryEnergy \|\| Tcut != lastTcut){
517	lastMaterial =aMaterial;
518	lastPrimaryEnergy = PrimEnergy;
519	lastTcut=Tcut;
520	listOfIndexOfAdjointEMModelInAction.clear();
521	listOfIsScatProjToProjCase.clear();
522	lastAdjointCSVsModelsAndElements.clear();
523	need_to_compute=true;
524
525	}
526	size_t ind=0;
527	if (!need_to_compute){
528	need_to_compute=true;
529	for (size_t i=0;i<listOfIndexOfAdjointEMModelInAction.size();i++){
530	size_t ind1=listOfIndexOfAdjointEMModelInAction[i];
531	if (aModel == listOfAdjointEMModel[ind1] && IsScatProjToProjCase == listOfIsScatProjToProjCase[i]){
532	need_to_compute=false;
533	CS_Vs_Element = lastAdjointCSVsModelsAndElements[ind];
534	}
535	ind++;
536	}
537	}
538
539	if (need_to_compute){
540	size_t ind_model=0;
541	for (size_t i=0;i<listOfAdjointEMModel.size();i++){
542	if (aModel == listOfAdjointEMModel[i]){
543	ind_model=i;
544	break;
545	}
546	}
547	G4double Tlow=Tcut;
548	if (!listOfAdjointEMModel[ind_model]->GetApplyCutInRange()) Tlow =listOfAdjointEMModel[ind_model]->GetLowEnergyLimit();
549	listOfIndexOfAdjointEMModelInAction.push_back(ind_model);
550	listOfIsScatProjToProjCase.push_back(IsScatProjToProjCase);
551	CS_Vs_Element.clear();
552	if (!aModel->GetUseMatrix()){
553	CS_Vs_Element.push_back(aModel->AdjointCrossSection(currentCouple,PrimEnergy,IsScatProjToProjCase));
554
555
556	}
557	else if (aModel->GetUseMatrixPerElement()){
558	size_t n_el = aMaterial->GetNumberOfElements();
559	if (aModel->GetUseOnlyOneMatrixForAllElements()){
560	G4AdjointCSMatrix* theCSMatrix;
561	if (IsScatProjToProjCase){
562	theCSMatrix=theAdjointCSMatricesForScatProjToProj[ind_model][0];
563	}
564	else theCSMatrix=theAdjointCSMatricesForProdToProj[ind_model][0];
565	G4double CS =0.;
566	if (PrimEnergy > Tlow)
567	CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow);
568	G4double factor=0.;
569	for (size_t i=0;i<n_el;i++){ //this could be computed only once
570	//size_t ind_el = aMaterial->GetElement(i)->GetIndex();
571	factor+=aMaterial->GetElement(i)->GetZ()*aMaterial->GetVecNbOfAtomsPerVolume()[i];
572	}
573	CS *=factor;
574	CS_Vs_Element.push_back(CS);
575
576	}
577	else {
578	for (size_t i=0;i<n_el;i++){
579	size_t ind_el = aMaterial->GetElement(i)->GetIndex();
580	//G4cout<<aMaterial->GetName()<<G4endl;
581	G4AdjointCSMatrix* theCSMatrix;
582	if (IsScatProjToProjCase){
583	theCSMatrix=theAdjointCSMatricesForScatProjToProj[ind_model][ind_el];
584	}
585	else theCSMatrix=theAdjointCSMatricesForProdToProj[ind_model][ind_el];
586	G4double CS =0.;
587	if (PrimEnergy > Tlow)
588	CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow);
589	//G4cout<<CS<<G4endl;
590	CS_Vs_Element.push_back(CS*(aMaterial->GetVecNbOfAtomsPerVolume()[i]));
591	}
592	}
593
594	}
595	else {
596	size_t ind_mat = aMaterial->GetIndex();
597	G4AdjointCSMatrix* theCSMatrix;
598	if (IsScatProjToProjCase){
599	theCSMatrix=theAdjointCSMatricesForScatProjToProj[ind_model][ind_mat];
600	}
601	else theCSMatrix=theAdjointCSMatricesForProdToProj[ind_model][ind_mat];
602	G4double CS =0.;
603	if (PrimEnergy > Tlow)
604	CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow);
605	CS_Vs_Element.push_back(CS);
606
607
608	}
609	lastAdjointCSVsModelsAndElements.push_back(CS_Vs_Element);
610
611	}
612
613
614	G4double CS=0;
615	for (size_t i=0;i<CS_Vs_Element.size();i++){
616	CS+=CS_Vs_Element[i]; //We could put the progressive sum of the CS instead of the CS of an element itself
617
618	}
619	return CS;
620	}
621	///////////////////////////////////////////////////////
622	//
623	G4Element* G4AdjointCSManager::SampleElementFromCSMatrices(G4Material* aMaterial,
624	G4VEmAdjointModel* aModel,
625	G4double PrimEnergy,
626	G4double Tcut,
627	G4bool IsScatProjToProjCase)
628	{ std::vector<G4double> CS_Vs_Element;
629	G4double CS = ComputeAdjointCS(aMaterial,aModel,PrimEnergy,Tcut,IsScatProjToProjCase,CS_Vs_Element);
630	G4double rand_var= G4UniformRand();
631	G4double SumCS=0.;
632	size_t ind=0;
633	for (size_t i=0;i<CS_Vs_Element.size();i++){
634	SumCS+=CS_Vs_Element[i];
635	if (rand_var<=SumCS/CS){
636	ind=i;
637	break;
638	}
639	}
640
641	return const_cast<G4Element*>(aMaterial->GetElement(ind));
642
643
644
645	}
646	///////////////////////////////////////////////////////
647	//
648	G4double G4AdjointCSManager::ComputeTotalAdjointCS(const G4MaterialCutsCouple* aCouple,
649	G4ParticleDefinition* aPartDef,
650	G4double Ekin)
651	{
652	G4double TotalCS=0.;
653
654	DefineCurrentMaterial(aCouple);
655
656
657	std::vector<G4double> CS_Vs_Element;
658	for (size_t i=0; i<listOfAdjointEMModel.size();i++){
659
660	G4double Tlow=0;
661	if (!listOfAdjointEMModel[i]->GetApplyCutInRange()) Tlow =listOfAdjointEMModel[i]->GetLowEnergyLimit();
662	else {
663	G4ParticleDefinition* theDirSecondPartDef =
664	GetForwardParticleEquivalent(listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition());
665	size_t idx=56;
666	if (theDirSecondPartDef->GetParticleName() == "gamma") idx = 0;
667	else if (theDirSecondPartDef->GetParticleName() == "e-") idx = 1;
668	else if (theDirSecondPartDef->GetParticleName() == "e+") idx = 2;
669	if (idx <56) {
670	const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
671	Tlow =(*aVec)[aCouple->GetIndex()];
672	}
673
674
675	}
676
677	if ( Ekin<=listOfAdjointEMModel[i]->GetHighEnergyLimit() && Ekin>=listOfAdjointEMModel[i]->GetLowEnergyLimit()){
678	if (aPartDef == listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectPrimaryParticleDefinition()){
679	TotalCS += ComputeAdjointCS(currentMaterial,
680	listOfAdjointEMModel[i],
681	Ekin, Tlow,true,CS_Vs_Element);
682	}
683	if (aPartDef == listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition()){
684	TotalCS += ComputeAdjointCS(currentMaterial,
685	listOfAdjointEMModel[i],
686	Ekin, Tlow,false, CS_Vs_Element);
687	}
688
689	}
690	}
691	return TotalCS;
692
693
694	}
695	///////////////////////////////////////////////////////
696	//
697	std::vector<G4AdjointCSMatrix*>
698	G4AdjointCSManager::BuildCrossSectionsMatricesForAGivenModelAndElement(G4VEmAdjointModel* aModel,G4int Z,G4int A,
699	G4int nbin_pro_decade)
700	{
701	G4AdjointCSMatrix* theCSMatForProdToProjBackwardScattering = new G4AdjointCSMatrix(false);
702	G4AdjointCSMatrix* theCSMatForScatProjToProjBackwardScattering = new G4AdjointCSMatrix(true);
703
704
705	//make the vector of primary energy of the adjoint particle, could try to make this just once ?
706
707	G4double EkinMin =aModel->GetLowEnergyLimit();
708	G4double EkinMaxForScat =aModel->GetHighEnergyLimit()*0.999;
709	G4double EkinMaxForProd =aModel->GetHighEnergyLimit()*0.999;
710	if (aModel->GetSecondPartOfSameType() )EkinMaxForProd =EkinMaxForProd/2.;
711
712
713	//Product to projectile backward scattering
714	//-----------------------------------------
715	G4double fE=std::pow(10.,1./nbin_pro_decade);
716	G4double E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
717	G4double E1=EkinMin;
718	while (E1 <EkinMaxForProd){
719	E1=std::max(EkinMin,E2);
720	E1=std::min(EkinMaxForProd,E1);
721	std::vector< std::vector< double>* > aMat= aModel->ComputeAdjointCrossSectionVectorPerAtomForSecond(E1,Z,A,nbin_pro_decade);
722	if (aMat.size()>=2) {
723	std::vector< double>* log_ESecVec=aMat[0];
724	std::vector< double>* log_CSVec=aMat[1];
725	G4double log_adjointCS=log_CSVec->back();
726	//normalise CSVec such that it becomes a probability vector
727	for (size_t j=0;j<log_CSVec->size();j++) {
728	if (j==0) (*log_CSVec)[j] = 0.;
729	else (log_CSVec)[j]=std::log(1.-std::exp((log_CSVec)[j]-log_adjointCS) +1e-50);
730	}
731	(log_CSVec)[log_CSVec->size()-1]=(log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
732	theCSMatForProdToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
733	}
734	E1=E2;
735	E2*=fE;
736	}
737
738	//Scattered projectile to projectile backward scattering
739	//-----------------------------------------
740
741	E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
742	E1=EkinMin;
743	while (E1 <EkinMaxForScat){
744	E1=std::max(EkinMin,E2);
745	E1=std::min(EkinMaxForScat,E1);
746	std::vector< std::vector< double>* > aMat= aModel->ComputeAdjointCrossSectionVectorPerAtomForScatProj(E1,Z,A,nbin_pro_decade);
747	if (aMat.size()>=2) {
748	std::vector< double>* log_ESecVec=aMat[0];
749	std::vector< double>* log_CSVec=aMat[1];
750	G4double log_adjointCS=log_CSVec->back();
751	//normalise CSVec such that it becomes a probability vector
752	for (size_t j=0;j<log_CSVec->size();j++) {
753	if (j==0) (*log_CSVec)[j] = 0.;
754	else (log_CSVec)[j]=std::log(1.-std::exp((log_CSVec)[j]-log_adjointCS)+1e-50);
755	}
756	(log_CSVec)[log_CSVec->size()-1]=(log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
757	theCSMatForScatProjToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
758	}
759	E1=E2;
760	E2*=fE;
761	}
762
763
764	std::vector<G4AdjointCSMatrix*> res;
765	res.clear();
766	res.push_back(theCSMatForProdToProjBackwardScattering);
767	res.push_back(theCSMatForScatProjToProjBackwardScattering);
768
769
770	/*
771	G4String file_name;
772	std::stringstream astream;
773	G4String str_Z;
774	astream<<Z;
775	astream>>str_Z;
776	theCSMatForProdToProjBackwardScattering->Write(aModel->GetName()+G4String("_CSMat_Z")+str_Z+"_ProdToProj.txt");
777	theCSMatForScatProjToProjBackwardScattering->Write(aModel->GetName()+G4String("_CSMat_Z")+str_Z+"_ScatProjToProj.txt");
778
779	*/
780
781
782	return res;
783
784
785	}
786	///////////////////////////////////////////////////////
787	//
788	std::vector<G4AdjointCSMatrix*>
789	G4AdjointCSManager::BuildCrossSectionsMatricesForAGivenModelAndMaterial(G4VEmAdjointModel* aModel,
790	G4Material* aMaterial,
791	G4int nbin_pro_decade)
792	{
793	G4AdjointCSMatrix* theCSMatForProdToProjBackwardScattering = new G4AdjointCSMatrix(false);
794	G4AdjointCSMatrix* theCSMatForScatProjToProjBackwardScattering = new G4AdjointCSMatrix(true);
795
796
797	//make the vector of primary energy of the adjoint particle, could try to make this just once ?
798
799	G4double EkinMin =aModel->GetLowEnergyLimit();
800	G4double EkinMaxForScat =aModel->GetHighEnergyLimit()*0.999;
801	G4double EkinMaxForProd =aModel->GetHighEnergyLimit()*0.999;
802	if (aModel->GetSecondPartOfSameType() )EkinMaxForProd =EkinMaxForProd/2.;
803
804
805
806
807
808
809
810	//Product to projectile backward scattering
811	//-----------------------------------------
812	G4double fE=std::pow(10.,1./nbin_pro_decade);
813	G4double E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
814	G4double E1=EkinMin;
815	while (E1 <EkinMaxForProd){
816	E1=std::max(EkinMin,E2);
817	E1=std::min(EkinMaxForProd,E1);
818	std::vector< std::vector< double>* > aMat= aModel->ComputeAdjointCrossSectionVectorPerVolumeForSecond(aMaterial,E1,nbin_pro_decade);
819	if (aMat.size()>=2) {
820	std::vector< double>* log_ESecVec=aMat[0];
821	std::vector< double>* log_CSVec=aMat[1];
822	G4double log_adjointCS=log_CSVec->back();
823
824	//normalise CSVec such that it becomes a probability vector
825	for (size_t j=0;j<log_CSVec->size();j++) {
826	//G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<G4endl;
827	if (j==0) (*log_CSVec)[j] = 0.;
828	else (log_CSVec)[j]=std::log(1.-std::exp((log_CSVec)[j]-log_adjointCS));
829	//G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<G4endl;
830	}
831	(log_CSVec)[log_CSVec->size()-1]=(log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
832	theCSMatForProdToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
833	}
834
835
836
837	E1=E2;
838	E2*=fE;
839	}
840
841	//Scattered projectile to projectile backward scattering
842	//-----------------------------------------
843
844	E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
845	E1=EkinMin;
846	while (E1 <EkinMaxForScat){
847	E1=std::max(EkinMin,E2);
848	E1=std::min(EkinMaxForScat,E1);
849	std::vector< std::vector< double>* > aMat= aModel->ComputeAdjointCrossSectionVectorPerVolumeForScatProj(aMaterial,E1,nbin_pro_decade);
850	if (aMat.size()>=2) {
851	std::vector< double>* log_ESecVec=aMat[0];
852	std::vector< double>* log_CSVec=aMat[1];
853	G4double log_adjointCS=log_CSVec->back();
854
855	for (size_t j=0;j<log_CSVec->size();j++) {
856	//G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<G4endl;
857	if (j==0) (*log_CSVec)[j] = 0.;
858	else (log_CSVec)[j]=std::log(1.-std::exp((log_CSVec)[j]-log_adjointCS));
859	//G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<G4endl;if (theAdjPartDef->GetParticleName() == "adj_gamma") return G4Gamma::Gamma();
860
861	}
862	(log_CSVec)[log_CSVec->size()-1]=(log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
863
864	theCSMatForScatProjToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
865	}
866	E1=E2;
867	E2*=fE;
868	}
869
870
871
872
873
874
875
876	std::vector<G4AdjointCSMatrix*> res;
877	res.clear();
878
879	res.push_back(theCSMatForProdToProjBackwardScattering);
880	res.push_back(theCSMatForScatProjToProjBackwardScattering);
881
882	/*
883	theCSMatForProdToProjBackwardScattering->Write(aModel->GetName()+"_CSMat_"+aMaterial->GetName()+"_ProdToProj.txt");
884	theCSMatForScatProjToProjBackwardScattering->Write(aModel->GetName()+"_CSMat_"+aMaterial->GetName()+"_ScatProjToProj.txt");
885	*/
886
887
888	return res;
889
890
891	}
892
893	///////////////////////////////////////////////////////
894	//
895	G4ParticleDefinition* G4AdjointCSManager::GetAdjointParticleEquivalent(G4ParticleDefinition* theFwdPartDef)
896	{
897	if (theFwdPartDef->GetParticleName() == "e-") return G4AdjointElectron::AdjointElectron();
898	else if (theFwdPartDef->GetParticleName() == "gamma") return G4AdjointGamma::AdjointGamma();
899	else if (theFwdPartDef->GetParticleName() == "proton") return G4AdjointProton::AdjointProton();
900	else if (theFwdPartDef ==theFwdIon) return theAdjIon;
901
902	return 0;
903	}
904	///////////////////////////////////////////////////////
905	//
906	G4ParticleDefinition* G4AdjointCSManager::GetForwardParticleEquivalent(G4ParticleDefinition* theAdjPartDef)
907	{
908	if (theAdjPartDef->GetParticleName() == "adj_e-") return G4Electron::Electron();
909	else if (theAdjPartDef->GetParticleName() == "adj_gamma") return G4Gamma::Gamma();
910	else if (theAdjPartDef->GetParticleName() == "adj_proton") return G4Proton::Proton();
911	else if (theAdjPartDef == theAdjIon) return theFwdIon;
912	return 0;
913	}
914	///////////////////////////////////////////////////////
915	//
916	void G4AdjointCSManager::DefineCurrentMaterial(const G4MaterialCutsCouple* couple)
917	{
918	if(couple != currentCouple) {
919	currentCouple = const_cast<G4MaterialCutsCouple*> (couple);
920	currentMaterial = const_cast<G4Material*> (couple->GetMaterial());
921	currentMatIndex = couple->GetIndex();
922	lastPartDefForCS =0;
923	LastEkinForCS =0;
924	LastCSCorrectionFactor =1.;
925	}
926	}
927
928	///////////////////////////////////////////////////////
929	//
930	void G4AdjointCSManager::DefineCurrentParticle(const G4ParticleDefinition* aPartDef)
931	{
932	if(aPartDef != currentParticleDef) {
933
934	currentParticleDef= const_cast< G4ParticleDefinition* > (aPartDef);
935	massRatio=1;
936	if (aPartDef == theAdjIon) massRatio = proton_mass_c2/aPartDef->GetPDGMass();
937	currentParticleIndex=1000000;
938	for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
939	if (aPartDef == theListOfAdjointParticlesInAction[i]) currentParticleIndex=i;
940	}
941
942	}
943	}
944
945
946
947	/////////////////////////////////////////////////////////////////////////////////////////////////
948	//
949	G4double G4AdjointCSManager::ComputeAdjointCS(G4double aPrimEnergy,G4AdjointCSMatrix*
950	anAdjointCSMatrix,G4double Tcut)
951	{
952	std::vector< double> *theLogPrimEnergyVector = anAdjointCSMatrix->GetLogPrimEnergyVector();
953	if (theLogPrimEnergyVector->size() ==0){
954	G4cout<<"No data are contained in the given AdjointCSMatrix!"<<G4endl;
955	G4cout<<"The s"<<G4endl;
956	return 0.;
957
958	}
959	G4double log_Tcut = std::log(Tcut);
960	G4double log_E =std::log(aPrimEnergy);
961
962	if (aPrimEnergy <= Tcut \|\| log_E > theLogPrimEnergyVector->back()) return 0.;
963
964
965
966	G4AdjointInterpolator* theInterpolator=G4AdjointInterpolator::GetInstance();
967
968	size_t ind =theInterpolator->FindPositionForLogVector(log_E,*theLogPrimEnergyVector);
969	G4double aLogPrimEnergy1,aLogPrimEnergy2;
970	G4double aLogCS1,aLogCS2;
971	G4double log01,log02;
972	std::vector< double>* aLogSecondEnergyVector1 =0;
973	std::vector< double>* aLogSecondEnergyVector2 =0;
974	std::vector< double>* aLogProbVector1=0;
975	std::vector< double>* aLogProbVector2=0;
976	std::vector< size_t>* aLogProbVectorIndex1=0;
977	std::vector< size_t>* aLogProbVectorIndex2=0;
978
979
980	anAdjointCSMatrix->GetData(ind, aLogPrimEnergy1,aLogCS1,log01, aLogSecondEnergyVector1,aLogProbVector1,aLogProbVectorIndex1);
981	anAdjointCSMatrix->GetData(ind+1, aLogPrimEnergy2,aLogCS2,log02, aLogSecondEnergyVector2,aLogProbVector2,aLogProbVectorIndex2);
982	if (anAdjointCSMatrix->IsScatProjToProjCase()){ //case where the Tcut plays a role
983	G4double log_minimum_prob1, log_minimum_prob2;
984	log_minimum_prob1=theInterpolator->InterpolateForLogVector(log_Tcut,aLogSecondEnergyVector1,aLogProbVector1);
985	log_minimum_prob2=theInterpolator->InterpolateForLogVector(log_Tcut,aLogSecondEnergyVector2,aLogProbVector2);
986	aLogCS1+= log_minimum_prob1;
987	aLogCS2+= log_minimum_prob2;
988	}
989
990	G4double log_adjointCS = theInterpolator->LinearInterpolation(log_E,aLogPrimEnergy1,aLogPrimEnergy2,aLogCS1,aLogCS2);
991	return std::exp(log_adjointCS);
992
993
994	}

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

Download in other formats: