Context Navigation

source: trunk/source/processes/electromagnetic/utils/src/G4EmModelManager.cc @ 1129

Last change on this file since 1129 was 1055, checked in by garnier, 15 years ago
maj sur la beta de geant 4.9.3
File size: 24.3 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: G4EmModelManager.cc,v 1.49 2009/04/17 10:35:32 vnivanch Exp $
27	// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
28	//
29	// -------------------------------------------------------------------
30	//
31	// GEANT4 Class file
32	//
33	//
34	// File name: G4EmModelManager
35	//
36	// Author: Vladimir Ivanchenko
37	//
38	// Creation date: 07.05.2002
39	//
40	// Modifications:
41	//
42	// 23-12-02 V.Ivanchenko change interface in order to move
43	// to cut per region
44	// 20-01-03 Migrade to cut per region (V.Ivanchenko)
45	// 24-01-03 Make models region aware (V.Ivanchenko)
46	// 13-02-03 The set of models is defined for region (V.Ivanchenko)
47	// 06-03-03 Fix in energy intervals for models (V.Ivanchenko)
48	// 13-04-03 Add startFromNull (V.Ivanchenko)
49	// 13-05-03 Add calculation of precise range (V.Ivanchenko)
50	// 16-07-03 Replace G4Material by G4MaterialCutCouple in dE/dx and CrossSection
51	// calculation (V.Ivanchenko)
52	// 21-07-03 Add UpdateEmModel method (V.Ivanchenko)
53	// 03-11-03 Substitute STL vector for G4RegionModels (V.Ivanchenko)
54	// 26-01-04 Fix in energy range conditions (V.Ivanchenko)
55	// 24-03-05 Remove check or IsInCharge (V.Ivanchenko)
56	// 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
57	// 18-08-05 Fix cut for e+e- pair production (V.Ivanchenko)
58	// 29-11-05 Add protection for arithmetic operations with cut=DBL_MAX (V.Ivanchenko)
59	// 20-01-06 Introduce G4EmTableType and reducing number of methods (VI)
60	// 13-05-06 Add GetModel by index method (VI)
61	// 15-03-07 Add maxCutInRange (V.Ivanchenko)
62	// 12-04-07 Add verbosity at destruction (V.Ivanchenko)
63	// 08-04-08 Fixed and simplified initialisation of G4RegionModel (VI)
64	//
65	// Class Description:
66	//
67	// It is the unified energy loss process it calculates the continuous
68	// energy loss for charged particles using a set of Energy Loss
69	// models valid for different energy regions. There are a possibility
70	// to create and access to dE/dx and range tables, or to calculate
71	// that information on fly.
72	// -------------------------------------------------------------------
73	//
74	//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75	//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76
77	#include "G4EmModelManager.hh"
78
79	//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
80
81	G4RegionModels::G4RegionModels(G4int nMod, std::vector<G4int>& indx,
82	G4DataVector& lowE, const G4Region* reg)
83	{
84	nModelsForRegion = nMod;
85	theListOfModelIndexes = new G4int [nModelsForRegion];
86	lowKineticEnergy = new G4double [nModelsForRegion+1];
87	for (G4int i=0; i<nModelsForRegion; i++) {
88	theListOfModelIndexes[i] = indx[i];
89	lowKineticEnergy[i] = lowE[i];
90	}
91	lowKineticEnergy[nModelsForRegion] = lowE[nModelsForRegion];
92	theRegion = reg;
93	}
94
95	//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
96
97	G4RegionModels::~G4RegionModels()
98	{
99	delete [] theListOfModelIndexes;
100	delete [] lowKineticEnergy;
101	}
102
103	//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
104
105	#include "G4Step.hh"
106	#include "G4ParticleDefinition.hh"
107	#include "G4PhysicsVector.hh"
108	#include "G4Gamma.hh"
109	#include "G4Positron.hh"
110	#include "G4MaterialCutsCouple.hh"
111	#include "G4ProductionCutsTable.hh"
112	#include "G4RegionStore.hh"
113	#include "G4Gamma.hh"
114	#include "G4Positron.hh"
115	#include "G4UnitsTable.hh"
116
117	//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
118
119	G4EmModelManager::G4EmModelManager():
120	nEmModels(0),
121	nRegions(0),
122	nCouples(0),
123	idxOfRegionModels(0),
124	setOfRegionModels(0),
125	minSubRange(0.1),
126	particle(0),
127	verboseLevel(0)
128	{
129	models.clear();
130	flucModels.clear();
131	regions.clear();
132	orderOfModels.clear();
133	maxCutInRange = 12.*cm;
134	maxSubCutInRange = 0.7*mm;
135	theGamma = G4Gamma::Gamma();
136	thePositron = G4Positron::Positron();
137	models.reserve(4);
138	flucModels.reserve(4);
139	regions.reserve(4);
140	orderOfModels.reserve(4);
141	isUsed.reserve(4);
142	}
143
144	//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
145
146	G4EmModelManager::~G4EmModelManager()
147	{
148	Clear();
149	}
150
151	//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
152
153	void G4EmModelManager::Clear()
154	{
155	if(1 < verboseLevel) {
156	G4cout << "G4EmModelManager::Clear()" << G4endl;
157	}
158
159	theCuts.clear();
160	theSubCuts.clear();
161	if(idxOfRegionModels) delete [] idxOfRegionModels;
162	if(setOfRegionModels && nRegions) {
163	for(G4int i=0; i<nRegions; i++) {
164	delete (setOfRegionModels[i]);
165	}
166	delete [] setOfRegionModels;
167	}
168	idxOfRegionModels = 0;
169	setOfRegionModels = 0;
170	}
171
172	//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
173
174	void G4EmModelManager::AddEmModel(G4int num, G4VEmModel* p,
175	G4VEmFluctuationModel* fm, const G4Region* r)
176	{
177	if(!p) {
178	G4cout << "G4EmModelManager::AddEmModel WARNING: no model defined."
179	<< G4endl;
180	return;
181	}
182	models.push_back(p);
183	flucModels.push_back(fm);
184	regions.push_back(r);
185	orderOfModels.push_back(num);
186	isUsed.push_back(0);
187	p->DefineForRegion(r);
188	nEmModels++;
189	}
190
191	//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
192
193	void G4EmModelManager::UpdateEmModel(const G4String& nam,
194	G4double emin, G4double emax)
195	{
196	if (nEmModels > 0) {
197	for(G4int i=0; i<nEmModels; i++) {
198	if(nam == models[i]->GetName()) {
199	models[i]->SetLowEnergyLimit(emin);
200	models[i]->SetHighEnergyLimit(emax);
201	break;
202	}
203	}
204	}
205	G4cout << "G4EmModelManager::UpdateEmModel WARNING: no model <"
206	<< nam << "> is found out"
207	<< G4endl;
208	}
209
210	//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
211
212	G4VEmModel* G4EmModelManager::GetModel(G4int i, G4bool ver)
213	{
214	G4VEmModel* m = 0;
215	if(i >= 0 && i < nEmModels) {m = models[i];}
216	else if(verboseLevel > 0 && ver) {
217	G4cout << "G4EmModelManager::GetModel WARNING: "
218	<< "index " << i << " is wrong Nmodels= "
219	<< nEmModels;
220	if(particle) G4cout << " for " << particle->GetParticleName();
221	G4cout<< G4endl;
222	}
223	return m;
224	}
225
226	//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
227
228	const G4DataVector* G4EmModelManager::Initialise(const G4ParticleDefinition* p,
229	const G4ParticleDefinition* sp,
230	G4double theMinSubRange,
231	G4int val)
232	{
233	verboseLevel = val;
234	G4String partname = p->GetParticleName();
235	if(1 < verboseLevel) {
236	G4cout << "G4EmModelManager::Initialise() for "
237	<< partname << G4endl;
238	}
239	// Are models defined?
240	if(!nEmModels) {
241	G4Exception("G4EmModelManager::Initialise without any model defined for "+partname);
242	}
243	particle = p;
244	secondaryParticle = sp;
245	minSubRange = theMinSubRange;
246
247	Clear();
248	G4RegionStore* regionStore = G4RegionStore::GetInstance();
249	const G4Region* world =
250	regionStore->GetRegion("DefaultRegionForTheWorld", false);
251
252	// Identify the list of regions with different set of models
253	nRegions = 1;
254	std::vector<const G4Region*> setr;
255	setr.push_back(world);
256	G4bool isWorld = false;
257
258	for (G4int ii=0; ii<nEmModels; ii++) {
259	const G4Region* r = regions[ii];
260	if ( r == 0 \|\| r == world) {
261	isWorld = true;
262	regions[ii] = world;
263	} else {
264	G4bool newRegion = true;
265	if (nRegions>1) {
266	for (G4int j=1; j<nRegions; j++) {
267	if ( r == setr[j] ) newRegion = false;
268	}
269	}
270	if (newRegion) {
271	setr.push_back(r);
272	nRegions++;
273	}
274	}
275	}
276
277	G4ProductionCutsTable* theCoupleTable=
278	G4ProductionCutsTable::GetProductionCutsTable();
279	G4int numOfCouples = theCoupleTable->GetTableSize();
280	if(nRegions > 1) idxOfRegionModels = new G4int[numOfCouples];
281	setOfRegionModels = new G4RegionModels*[nRegions];
282
283	std::vector<G4int> modelAtRegion(nEmModels);
284	std::vector<G4int> modelOrd(nEmModels);
285	G4DataVector eLow(nEmModels+1);
286	G4DataVector eHigh(nEmModels);
287
288	// Order models for regions
289	for (G4int reg=0; reg<nRegions; reg++) {
290	const G4Region* region = setr[reg];
291	G4int n = 0;
292
293	if(isWorld \|\| 0 < reg) {
294
295	for (G4int ii=0; ii<nEmModels; ii++) {
296
297	G4VEmModel* model = models[ii];
298	if ( region == regions[ii] ) {
299
300	G4double tmin = model->LowEnergyLimit();
301	G4double tmax = model->HighEnergyLimit();
302	G4int ord = orderOfModels[ii];
303	G4bool push = true;
304	G4bool insert = false;
305	G4int idx = n;
306
307	if(1 < verboseLevel) {
308	G4cout << "Model #" << ii
309	<< " <" << model->GetName() << "> for region <";
310	if (region) G4cout << region->GetName();
311	G4cout << "> "
312	<< " tmin(MeV)= " << tmin/MeV
313	<< "; tmax(MeV)= " << tmax/MeV
314	<< "; order= " << ord
315	<< G4endl;
316	}
317
318	if(n > 0) {
319
320	// extend energy range to previous models
321	tmin = std::min(tmin, eHigh[n-1]);
322	tmax = std::max(tmax, eLow[0]);
323	//G4cout << "tmin= " << tmin << " tmax= "
324	// << tmax << " ord= " << ord <<G4endl;
325	// empty energy range
326	if( tmax - tmin <= eV) push = false;
327	// low-energy model
328	else if (tmax == eLow[0]) {
329	push = false;
330	insert = true;
331	idx = 0;
332	// resolve intersections
333	} else if(tmin < eHigh[n-1]) {
334	// compare order
335	for(G4int k=0; k<n; k++) {
336	// new model has lower application
337	if(ord >= modelOrd[k]) {
338	if(tmin < eHigh[k] && tmin >= eLow[k]) tmin = eHigh[k];
339	if(tmax <= eHigh[k] && tmax > eLow[k]) tmax = eLow[k];
340	if(tmax > eHigh[k] && tmin < eLow[k]) {
341	if(tmax - eHigh[k] > eLow[k] - tmin) tmin = eHigh[k];
342	else tmax = eLow[k];
343	}
344	if( tmax - tmin <= eV) {
345	push = false;
346	break;
347	}
348	}
349	}
350	//G4cout << "tmin= " << tmin << " tmax= "
351	// << tmax << " push= " << push << " idx= " << idx <<G4endl;
352	if(push) {
353	if (tmax == eLow[0]) {
354	push = false;
355	insert = true;
356	idx = 0;
357	// continue resolve intersections
358	} else if(tmin < eHigh[n-1]) {
359	// last energy interval
360	if(tmin > eLow[n-1] && tmax >= eHigh[n-1]) {
361	eHigh[n-1] = tmin;
362	// first energy interval
363	} else if(tmin <= eLow[0] && tmax < eHigh[0]) {
364	eLow[0] = tmax;
365	push = false;
366	insert = true;
367	idx = 0;
368	} else {
369	// find energy interval to replace
370	for(G4int k=0; k<n; k++) {
371	if(tmin <= eLow[k] && tmax >= eHigh[k]) {
372	push = false;
373	modelAtRegion[k] = ii;
374	modelOrd[k] = ord;
375	isUsed[ii] = 1;
376	}
377	}
378	}
379	}
380	}
381	}
382	}
383	if(insert) {
384	for(G4int k=n-1; k>=idx; k--) {
385	modelAtRegion[k+1] = modelAtRegion[k];
386	modelOrd[k+1] = modelOrd[k];
387	eLow[k+1] = eLow[k];
388	eHigh[k+1] = eHigh[k];
389	}
390	}
391	//G4cout << "push= " << push << " insert= " << insert
392	//<< " idx= " << idx <<G4endl;
393	if (push \|\| insert) {
394	n++;
395	modelAtRegion[idx] = ii;
396	modelOrd[idx] = ord;
397	eLow[idx] = tmin;
398	eHigh[idx] = tmax;
399	isUsed[ii] = 1;
400	}
401	}
402	}
403	} else {
404	n = 1;
405	models.push_back(0);
406	modelAtRegion.push_back(nEmModels);
407	eLow.push_back(0.0);
408	eHigh.push_back(DBL_MAX);
409	}
410	eLow[0] = 0.0;
411	eLow[n] = eHigh[n-1];
412
413	if(1 < verboseLevel) {
414	G4cout << "New G4RegionModels set with " << n << " models for region <";
415	if (region) G4cout << region->GetName();
416	G4cout << "> Elow(MeV)= ";
417	for(G4int ii=0; ii<=n; ii++) {G4cout << eLow[ii]/MeV << " ";}
418	G4cout << G4endl;
419	}
420	G4RegionModels* rm = new G4RegionModels(n, modelAtRegion, eLow, region);
421	setOfRegionModels[reg] = rm;
422	}
423
424	currRegionModel = setOfRegionModels[0];
425
426	// Access to materials and build cuts
427	for(G4int i=0; i<numOfCouples; i++) {
428
429	const G4MaterialCutsCouple* couple =
430	theCoupleTable->GetMaterialCutsCouple(i);
431	const G4Material* material = couple->GetMaterial();
432	const G4ProductionCuts* pcuts = couple->GetProductionCuts();
433
434	G4int reg = 0;
435	if(nRegions > 1) {
436	reg = nRegions;
437	do {reg--;} while (reg>0 && pcuts != (setr[reg]->GetProductionCuts()));
438	idxOfRegionModels[i] = reg;
439	}
440	if(1 < verboseLevel) {
441	G4cout << "G4EmModelManager::Initialise() for "
442	<< material->GetName()
443	<< " indexOfCouple= " << i
444	<< " indexOfRegion= " << reg
445	<< G4endl;
446	}
447
448	G4double cut = DBL_MAX;
449	G4double subcut = DBL_MAX;
450	if(secondaryParticle) {
451	size_t idx = 1;
452	if( secondaryParticle == theGamma ) idx = 0;
453	cut = (*theCoupleTable->GetEnergyCutsVector(idx))[i];
454	if( secondaryParticle == thePositron && cut < DBL_MAX )
455	cut += (*theCoupleTable->GetEnergyCutsVector(2))[i] +
456	2.0*electron_mass_c2;
457
458	// compute subcut
459	if( cut < DBL_MAX ) subcut = minSubRange*cut;
460	if(pcuts->GetProductionCut(idx) < maxCutInRange) {
461	G4double tcutmax =
462	theCoupleTable->ConvertRangeToEnergy(secondaryParticle,
463	material,maxSubCutInRange);
464	if(tcutmax < subcut) subcut = tcutmax;
465	}
466	}
467
468	G4int nm = setOfRegionModels[reg]->NumberOfModels();
469	for(G4int j=0; j<nm; j++) {
470
471	G4VEmModel* model = models[setOfRegionModels[reg]->ModelIndex(j)];
472
473	G4double tcutmin = model->MinEnergyCut(particle, couple);
474
475	if(cut < tcutmin) cut = tcutmin;
476	if(subcut < tcutmin) subcut = tcutmin;
477	if(1 < verboseLevel) {
478	G4cout << "The model # " << j
479	<< "; tcutmin(MeV)= " << tcutmin/MeV
480	<< "; tcut(MeV)= " << cut/MeV
481	<< "; tsubcut(MeV)= " << subcut/MeV
482	<< " for " << particle->GetParticleName()
483	<< G4endl;
484	}
485	}
486	theCuts.push_back(cut);
487	theSubCuts.push_back(subcut);
488	}
489
490	for(G4int jj=0; jj<nEmModels; jj++) {
491	if(1 == isUsed[jj]) {
492	models[jj]->Initialise(particle, theCuts);
493	if(flucModels[jj]) flucModels[jj]->InitialiseMe(particle);
494	}
495	}
496
497	if(1 < verboseLevel) {
498	G4cout << "G4EmModelManager for " << particle->GetParticleName()
499	<< " is initialised; nRegions= " << nRegions
500	<< G4endl;
501	}
502
503	return &theCuts;
504	}
505
506	//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
507
508	void G4EmModelManager::FillDEDXVector(G4PhysicsVector* aVector,
509	const G4MaterialCutsCouple* couple,
510	G4EmTableType tType)
511	{
512
513	G4double e;
514
515	size_t i = couple->GetIndex();
516	G4double cut = theCuts[i];
517	G4double subcut = 0.0;
518
519	if(fTotal == tType) cut = DBL_MAX;
520	else if(fSubRestricted == tType) subcut = theSubCuts[i];
521
522	if(1 < verboseLevel) {
523	G4cout << "G4EmModelManager::FillDEDXVector() for "
524	<< couple->GetMaterial()->GetName()
525	<< " cut(MeV)= " << cut
526	<< " subcut(MeV)= " << subcut
527	<< " Type " << tType
528	<< " for " << particle->GetParticleName()
529	<< G4endl;
530	}
531
532	G4int reg = 0;
533	if(nRegions > 1) reg = idxOfRegionModels[i];
534	const G4RegionModels* regModels = setOfRegionModels[reg];
535	G4int nmod = regModels->NumberOfModels();
536
537	// vectors to provide continues dE/dx
538	G4DataVector factor(nmod);
539	G4DataVector eLow(nmod+1);
540	G4DataVector dedxLow(nmod);
541	G4DataVector dedxHigh(nmod);
542
543	if(1 < verboseLevel) {
544	G4cout << "There are " << nmod << " models for "
545	<< couple->GetMaterial()->GetName()
546	<< " at the region #" << reg
547	<< G4endl;
548	}
549
550	// calculate factors to provide continuity of energy loss
551
552
553	factor[0] = 1.0;
554	G4int j;
555
556	G4int totBinsLoss = aVector->GetVectorLength();
557
558	dedxLow[0] = 0.0;
559	eLow[0] = 0.0;
560
561	e = regModels->LowEdgeEnergy(1);
562	eLow[1] = e;
563	G4VEmModel* model = models[regModels->ModelIndex(0)];
564	dedxHigh[0] = 0.0;
565	if(model && cut > subcut) {
566	dedxHigh[0] = model->ComputeDEDX(couple,particle,e,cut);
567	if(subcut > 0.0) {
568	dedxHigh[0] -= model->ComputeDEDX(couple,particle,e,subcut);
569	}
570	}
571	if(nmod > 1) {
572	for(j=1; j<nmod; j++) {
573
574	e = regModels->LowEdgeEnergy(j);
575	eLow[j] = e;
576	G4int idx = regModels->ModelIndex(j);
577
578	dedxLow[j] = models[idx]->ComputeDEDX(couple,particle,e,cut);
579	if(subcut > 0.0) {
580	dedxLow[j] -= models[idx]->ComputeDEDX(couple,particle,e,subcut);
581	}
582	if(subcut == cut) dedxLow[j] = 0.0;
583
584	e = regModels->LowEdgeEnergy(j+1);
585	eLow[j+1] = e;
586	dedxHigh[j] = models[idx]->ComputeDEDX(couple,particle,e,cut);
587	if(subcut > 0.0) {
588	dedxHigh[j] -= models[idx]->ComputeDEDX(couple,particle,e,subcut);
589	}
590	if(subcut == cut) dedxHigh[j] = 0.0;
591	}
592	if(1 < verboseLevel) {
593	G4cout << " model #0"
594	<< " dedx(" << eLow[0] << ")= " << dedxLow[0]
595	<< " dedx(" << eLow[1] << ")= " << dedxHigh[0]
596	<< G4endl;
597	}
598
599	for(j=1; j<nmod; j++) {
600	if(dedxLow[j] > 0.0) {
601	factor[j] = (dedxHigh[j-1]/dedxLow[j] - 1.0)*eLow[j];
602	} else factor[j] = 0.0;
603	if(1 < verboseLevel) {
604	G4cout << " model #" << j
605	<< " dedx(" << eLow[j] << ")= " << dedxLow[j]
606	<< " dedx(" << eLow[j+1] << ")= " << dedxHigh[j]
607	<< " factor= " << factor[j]/eLow[j]
608	<< G4endl;
609	}
610	}
611
612	if(2 < verboseLevel) {
613	G4cout << "Loop over " << totBinsLoss << " bins start " << G4endl;
614	}
615	}
616
617	// Calculate energy losses vector
618	for(j=0; j<totBinsLoss; j++) {
619
620	G4double e = aVector->GetLowEdgeEnergy(j);
621	G4double fac = 1.0;
622
623	// Choose a model of energy losses
624	G4int k = 0;
625	if (nmod > 1 && e > eLow[1]) {
626	do {
627	k++;
628	fac *= (1.0 + factor[k]/e);
629	} while (k+1 < nmod && e > eLow[k+1]);
630	}
631
632	model = models[regModels->ModelIndex(k)];
633	G4double dedx = 0.0;
634	G4double dedx0 = 0.0;
635	if(model && cut > subcut) {
636	dedx = model->ComputeDEDX(couple,particle,e,cut);
637	dedx0 = dedx;
638	if(subcut > 0.0) dedx -= model->ComputeDEDX(couple,particle,e,subcut);
639	dedx *= fac;
640	}
641
642	if(dedx < 0.0) dedx = 0.0;
643	if(2 < verboseLevel) {
644	G4cout << "Material= " << couple->GetMaterial()->GetName()
645	<< " E(MeV)= " << e/MeV
646	<< " dEdx(MeV/mm)= " << dedx*mm/MeV
647	<< " dEdx0(MeV/mm)= " << dedx0*mm/MeV
648	<< " fac= " << fac
649	<< G4endl;
650	}
651	aVector->PutValue(j, dedx);
652	}
653	}
654
655	//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
656
657	void G4EmModelManager::FillLambdaVector(G4PhysicsVector* aVector,
658	const G4MaterialCutsCouple* couple,
659	G4bool startFromNull,
660	G4EmTableType tType)
661	{
662	G4double e;
663
664	size_t i = couple->GetIndex();
665	G4double cut = theCuts[i];
666	G4double tmax = DBL_MAX;
667	if (fSubRestricted == tType) {
668	tmax = cut;
669	cut = theSubCuts[i];
670	}
671
672	if(1 < verboseLevel) {
673	G4cout << "G4EmModelManager::FillLambdaVector() for particle "
674	<< particle->GetParticleName()
675	<< " in " << couple->GetMaterial()->GetName()
676	<< " Ecut(MeV)= " << cut
677	<< " Emax(MeV)= " << tmax
678	<< " Type " << tType
679	<< G4endl;
680	}
681
682	G4int reg = 0;
683	if(nRegions > 1) reg = idxOfRegionModels[i];
684	const G4RegionModels* regModels = setOfRegionModels[reg];
685	G4int nmod = regModels->NumberOfModels();
686
687	// vectors to provide continues dE/dx
688	G4DataVector factor(nmod);
689	G4DataVector eLow(nmod+1);
690	G4DataVector sigmaLow(nmod);
691	G4DataVector sigmaHigh(nmod);
692
693	if(2 < verboseLevel) {
694	G4cout << "There are " << nmod << " models for "
695	<< couple->GetMaterial()->GetName() << G4endl;
696	}
697
698	// calculate factors to provide continuity of energy loss
699	factor[0] = 1.0;
700	G4int j;
701	G4int totBinsLambda = aVector->GetVectorLength();
702
703	sigmaLow[0] = 0.0;
704	eLow[0] = 0.0;
705
706	e = regModels->LowEdgeEnergy(1);
707	eLow[1] = e;
708	G4VEmModel* model = models[regModels->ModelIndex(0)];
709	sigmaHigh[0] = 0.0;
710	if(model) sigmaHigh[0] = model->CrossSection(couple,particle,e,cut,tmax);
711
712	if(2 < verboseLevel) {
713	G4cout << "### For material " << couple->GetMaterial()->GetName()
714	<< " " << nmod
715	<< " models"
716	<< " Ecut(MeV)= " << cut/MeV
717	<< " Emax(MeV)= " << e/MeV
718	<< " nbins= " << totBinsLambda
719	<< G4endl;
720	G4cout << " model #0 eUp= " << e
721	<< " sigmaUp= " << sigmaHigh[0] << G4endl;
722	}
723
724
725	if(nmod > 1) {
726
727	for(j=1; j<nmod; j++) {
728
729	e = regModels->LowEdgeEnergy(j);
730	eLow[j] = e;
731	G4int idx = regModels->ModelIndex(j);
732
733	sigmaLow[j] = models[idx]->CrossSection(couple,particle,e,cut,tmax);
734	e = regModels->LowEdgeEnergy(j+1);
735	eLow[j+1] = e;
736	sigmaHigh[j] = models[idx]->CrossSection(couple,particle,e,cut,tmax);
737	}
738	if(1 < verboseLevel) {
739	G4cout << " model #0"
740	<< " sigma(" << eLow[0] << ")= " << sigmaLow[0]
741	<< " sigma(" << eLow[1] << ")= " << sigmaHigh[0]
742	<< G4endl;
743	}
744	for(j=1; j<nmod; j++) {
745	if(sigmaLow[j] > 0.0) {
746	factor[j] = (sigmaHigh[j-1]/sigmaLow[j] - 1.0)*eLow[j];
747	} else factor[j] = 0.0;
748	if(1 < verboseLevel) {
749	G4cout << " model #" << j
750	<< " sigma(" << eLow[j] << ")= " << sigmaLow[j]
751	<< " sigma(" << eLow[j+1] << ")= " << sigmaHigh[j]
752	<< " factor= " << factor[j]/eLow[j]
753	<< G4endl;
754	}
755	}
756	}
757
758	// Calculate lambda vector
759	for(j=0; j<totBinsLambda; j++) {
760
761	e = aVector->GetLowEdgeEnergy(j);
762
763	// Choose a model of energy losses
764	G4int k = 0;
765	G4double fac = 1.0;
766	if (nmod > 1 && e > eLow[1]) {
767	do {
768	k++;
769	fac *= (1.0 + factor[k]/e);
770	} while ( k+1 < nmod && e > eLow[k+1] );
771	}
772
773	model = models[regModels->ModelIndex(k)];
774	G4double cross = 0.0;
775	if(model) cross = model->CrossSection(couple,particle,e,cut,tmax)*fac;
776	if(j==0 && startFromNull) cross = 0.0;
777
778	if(2 < verboseLevel) {
779	G4cout << "FillLambdaVector: " << j << ". e(MeV)= " << e/MeV
780	<< " cross(1/mm)= " << cross*mm
781	<< " fac= " << fac << " k= " << k
782	<< " model= " << regModels->ModelIndex(k)
783	<< G4endl;
784	}
785	if(cross < 0.0) cross = 0.0;
786
787	aVector->PutValue(j, cross);
788	}
789	}
790
791	//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
792
793	void G4EmModelManager::DumpModelList(G4int verb)
794	{
795	if(verb == 0) return;
796	for(G4int i=0; i<nRegions; i++) {
797	G4RegionModels* r = setOfRegionModels[i];
798	const G4Region* reg = r->Region();
799	// if(verb > -1 \|\| nRegions > 1) {
800	// }
801	G4int n = r->NumberOfModels();
802	if(verb > -1 \|\| n > 0) {
803	G4cout << " ===== EM models for the G4Region " << reg->GetName()
804	<< " ======" << G4endl;;
805	for(G4int j=0; j<n; j++) {
806	const G4VEmModel* m = models[r->ModelIndex(j)];
807	G4cout << std::setw(20);
808	G4cout << m->GetName() << " : Emin= "
809	<< std::setw(10) << G4BestUnit(r->LowEdgeEnergy(j),"Energy")
810	<< " Emax= "
811	<< G4BestUnit(r->LowEdgeEnergy(j+1),"Energy")
812	<< G4endl;
813	}
814	}
815	}
816	}
817
818	//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

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

Download in other formats: