Context Navigation

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

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

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

Download in other formats: