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

Last change on this file since 1193 was 1055, checked in by garnier, 17 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
81G4RegionModels::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
97G4RegionModels::~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
119G4EmModelManager::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
146G4EmModelManager::~G4EmModelManager()
147{
148 Clear();
149}
150
151//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
152
153void 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
174void 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
193void 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
212G4VEmModel* 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
228const 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
508void 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
657void 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
793void 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.