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

Last change on this file since 991 was 961, checked in by garnier, 17 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
80G4RegionModels::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
96G4RegionModels::~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
118G4EmModelManager::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
140G4EmModelManager::~G4EmModelManager()
141{
142 Clear();
143}
144
145//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
146
147void 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
168void 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
186void 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
205G4VEmModel* 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
221const 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
512void 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
660void 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
795void 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.