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

Last change on this file since 953 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 21.7 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// $Id: G4EmModelManager.cc,v 1.40 2007/11/09 11:35:54 vnivanch Exp $
27// GEANT4 tag $Name: $
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>& list, G4DataVector& lowE)
81{
82 nModelsForRegion = nMod;
83 theListOfModelIndexes = new G4int [nModelsForRegion];
84 lowKineticEnergy = new G4double [nModelsForRegion];
85 for (G4int i=0; i<nModelsForRegion; i++) {
86 theListOfModelIndexes[i] = list[i];
87 lowKineticEnergy[i] = lowE[i];
88 }
89}
90
91//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92
93G4RegionModels::~G4RegionModels()
94{
95 delete [] theListOfModelIndexes;
96 delete [] lowKineticEnergy;
97}
98
99//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100
101#include "G4LossTableManager.hh"
102#include "G4Step.hh"
103#include "G4ParticleDefinition.hh"
104#include "G4PhysicsVector.hh"
105#include "G4Gamma.hh"
106#include "G4Positron.hh"
107#include "G4MaterialCutsCouple.hh"
108#include "G4ProductionCutsTable.hh"
109#include "G4Region.hh"
110#include "G4RegionStore.hh"
111#include "G4Gamma.hh"
112#include "G4Positron.hh"
113
114//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
115
116G4EmModelManager::G4EmModelManager():
117 nEmModels(0),
118 nRegions(0),
119 nCouples(0),
120 idxOfRegionModels(0),
121 setOfRegionModels(0),
122 minSubRange(0.1),
123 particle(0),
124 verboseLevel(0)
125{
126 models.clear();
127 flucModels.clear();
128 regions.clear();
129 orderOfModels.clear();
130 upperEkin.clear();
131 maxCutInRange = 12.*cm;
132 maxSubCutInRange = 0.7*mm;
133 theGamma = G4Gamma::Gamma();
134 thePositron = G4Positron::Positron();
135}
136
137//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
138
139G4EmModelManager::~G4EmModelManager()
140{
141 G4int i,j;
142 Clear();
143 if(1 < verboseLevel) {
144 G4cout << "G4EmModelManager:: delete models ";
145 if(particle) G4cout << " for " << particle->GetParticleName();
146 G4cout << " nModels=" << nEmModels <<G4endl;
147 }
148
149 for(i = 0; i<nEmModels; i++) {
150 orderOfModels[i] = 1;
151 }
152 for(i = 0; i<nEmModels; i++) {
153 if (orderOfModels[i]) {
154 orderOfModels[i] = 0;
155 for(j = i+1; j<nEmModels; j++) {
156 if(models[i] == models[j]) orderOfModels[j] = 0;
157 }
158 G4String nam = models[i]->GetName();
159 if(nam != "PAI" && nam != "PAIModel" ) delete models[i];
160 }
161 }
162 for(i = 0; i<nEmModels; i++) {
163 orderOfModels[i] = 1;
164 }
165 for(i = 0; i<nEmModels; i++) {
166 if (orderOfModels[i]) {
167 orderOfModels[i] = 0;
168 for(j = i+1; j<nEmModels; j++) {
169 if(flucModels[i] == flucModels[j]) orderOfModels[j] = 0;
170 }
171 delete flucModels[i];
172 }
173 }
174 if(1 < verboseLevel)
175 G4cout << "G4EmModelManager:: models are deleted!" << G4endl;
176}
177
178//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
179
180void G4EmModelManager::Clear()
181{
182 if(1 < verboseLevel) {
183 G4cout << "G4EmModelManager::Clear()";
184 if(particle) G4cout << " for " << particle->GetParticleName();
185 G4cout << G4endl;
186 }
187
188 theCuts.clear();
189 theSubCuts.clear();
190 upperEkin.clear();
191 if(idxOfRegionModels) delete [] idxOfRegionModels;
192 if(setOfRegionModels && nRegions) {
193 for(G4int i=0; i<nRegions; i++) {
194 delete (setOfRegionModels[i]);
195 }
196 delete [] setOfRegionModels;
197 }
198 idxOfRegionModels = 0;
199 setOfRegionModels = 0;
200}
201
202//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
203
204void G4EmModelManager::AddEmModel(G4int num, G4VEmModel* p,
205 G4VEmFluctuationModel* fm, const G4Region* r)
206{
207 if(!p) {
208 G4cout << "G4EmModelManager::AddEmModel WARNING: no model defined."
209 << G4endl;
210 return;
211 }
212 models.push_back(p);
213 flucModels.push_back(fm);
214 regions.push_back(r);
215 orderOfModels.push_back(num);
216 p->DefineForRegion(r);
217 if (nEmModels>0) {
218 G4int idx = nEmModels;
219 do {idx--;} while (idx && num < orderOfModels[idx]);
220 if (num >= orderOfModels[idx] && num <= orderOfModels[idx+1]) idx++;
221 if (idx < nEmModels) {
222 models[nEmModels] = models[idx];
223 flucModels[nEmModels] = flucModels[idx];
224 regions[nEmModels] = regions[idx];
225 orderOfModels[nEmModels] = orderOfModels[idx];
226 models[idx] = p;
227 flucModels[idx] = fm;
228 regions[idx] = r;
229 orderOfModels[idx] = num;
230 }
231 }
232 nEmModels++;
233}
234
235//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
236
237void G4EmModelManager::UpdateEmModel(const G4String& nam,
238 G4double emin, G4double emax)
239{
240 if (nEmModels) {
241 for(G4int i=0; i<nEmModels; i++) {
242 if(nam == models[i]->GetName()) {
243 models[i]->SetLowEnergyLimit(emin);
244 models[i]->SetHighEnergyLimit(emax);
245 break;
246 }
247 }
248 }
249 G4cout << "G4EmModelManager::UpdateEmModel WARNING: no model <"
250 << nam << "> is found out"
251 << G4endl;
252}
253
254//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
255
256G4VEmModel* G4EmModelManager::GetModel(G4int i)
257{
258 G4VEmModel* m = 0;
259 if(i >= 0 && i < nEmModels) m = models[i];
260 else if(verboseLevel > 0)
261 G4cout << "G4EmModelManager::GetModel WARNING: "
262 << "index " << i << " is wrong Nmodels= "
263 << nEmModels << G4endl;
264 return m;
265}
266
267//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
268
269const G4DataVector* G4EmModelManager::Initialise(const G4ParticleDefinition* p,
270 const G4ParticleDefinition* sp,
271 G4double theMinSubRange,
272 G4int val)
273{
274 verboseLevel = val;
275 if(1 < verboseLevel) {
276 G4cout << "G4EmModelManager::Initialise() for "
277 << p->GetParticleName()
278 << G4endl;
279 }
280 // Are models defined?
281 if(!nEmModels) {
282 G4Exception("G4EmModelManager::Initialise without any model defined");
283 }
284 particle = p;
285 secondaryParticle = sp;
286 minSubRange = theMinSubRange;
287
288 Clear();
289 G4RegionStore* regionStore = G4RegionStore::GetInstance();
290 const G4Region* world =
291 regionStore->GetRegion("DefaultRegionForTheWorld", false);
292
293 // Identify the list of regions with different set of models
294 nRegions = 1;
295 std::vector<const G4Region*> set;
296 set.push_back(world);
297 G4bool isWorld = false;
298
299 for (G4int ii=0; ii<nEmModels; ii++) {
300 const G4Region* r = regions[ii];
301 if ( r == 0 || r == world) {
302 isWorld = true;
303 regions[ii] = world;
304 } else {
305 G4bool newRegion = true;
306 if (nRegions>1) {
307 for (G4int j=1; j<nRegions; j++) {
308 if ( r == set[j] ) newRegion = false;
309 }
310 }
311 if (newRegion) {
312 set.push_back(r);
313 nRegions++;
314 }
315 }
316 }
317
318 G4ProductionCutsTable* theCoupleTable=
319 G4ProductionCutsTable::GetProductionCutsTable();
320 G4int numOfCouples = theCoupleTable->GetTableSize();
321 idxOfRegionModels = new G4int[numOfCouples+1];
322 idxOfRegionModels[numOfCouples] = 0;
323 setOfRegionModels = new G4RegionModels*[nRegions];
324 upperEkin.resize(nEmModels);
325
326 // Order models for regions
327 for (G4int reg=0; reg<nRegions; reg++) {
328 const G4Region* region = set[reg];
329
330 G4int n = 0;
331
332 std::vector<G4int> modelAtRegion;
333 G4DataVector eLow;
334 G4DataVector eHigh;
335 modelAtRegion.clear();
336 eLow.clear();
337 eHigh.clear();
338
339 if(isWorld || 0 < reg) {
340
341 for (G4int ii=0; ii<nEmModels; ii++) {
342
343 G4VEmModel* model = models[ii];
344 if ( region == regions[ii] ) {
345
346 G4double tmin = model->LowEnergyLimit();
347 G4double tmax = model->HighEnergyLimit();
348 if (n>0) tmin = std::max(tmin, eHigh[n-1]);
349
350 if(1 < verboseLevel) {
351 G4cout << "Model #" << ii
352 << " <" << model->GetName() << "> for region <";
353 if (region) G4cout << region->GetName();
354 G4cout << "> "
355 << " tmin(MeV)= " << tmin/MeV
356 << "; tmax(MeV)= " << tmax/MeV
357 << "; order= " << orderOfModels[ii]
358 << G4endl;
359 }
360
361 if (tmin < tmax) {
362 modelAtRegion.push_back(ii);
363 eLow.push_back(tmin);
364 eHigh.push_back(tmax);
365 upperEkin[ii] = tmax;
366 n++;
367 }
368 }
369 }
370 } else {
371 n = 1;
372 models.push_back(0);
373 modelAtRegion.push_back(nEmModels);
374 eLow.push_back(0.0);
375 eHigh.push_back(DBL_MAX);
376 upperEkin.push_back(DBL_MAX);
377 }
378 eLow[0] = 0.0;
379
380 if(1 < verboseLevel) {
381 G4cout << "New G4RegionModels set with " << n << " models for region <";
382 if (region) G4cout << region->GetName();
383 G4cout << "> Elow(MeV)= ";
384 for(G4int ii=0; ii<n; ii++) {G4cout << eLow[ii]/MeV << " ";}
385 G4cout << G4endl;
386 }
387 G4RegionModels* rm = new G4RegionModels(n, modelAtRegion, eLow);
388 setOfRegionModels[reg] = rm;
389 }
390
391 // Access to materials and build cuts
392
393 for(G4int i=0; i<numOfCouples; i++) {
394
395 const G4MaterialCutsCouple* couple =
396 theCoupleTable->GetMaterialCutsCouple(i);
397 const G4Material* material = couple->GetMaterial();
398 const G4ProductionCuts* pcuts = couple->GetProductionCuts();
399
400 G4int reg = nRegions;
401 do {reg--;} while (reg>0 && pcuts != (set[reg]->GetProductionCuts()));
402 idxOfRegionModels[i] = reg;
403
404 if(1 < verboseLevel) {
405 G4cout << "G4EmModelManager::Initialise() for "
406 << material->GetName()
407 << " indexOfCouple= " << i
408 << " indexOfRegion= " << reg
409 << G4endl;
410 }
411
412 G4double cut = DBL_MAX;
413 G4double subcut = DBL_MAX;
414 if(secondaryParticle) {
415 size_t idx = 1;
416 if( secondaryParticle == theGamma ) idx = 0;
417 cut = (*theCoupleTable->GetEnergyCutsVector(idx))[i];
418 if( secondaryParticle == thePositron && cut < DBL_MAX )
419 cut += (*theCoupleTable->GetEnergyCutsVector(2))[i] +
420 2.0*electron_mass_c2;
421
422 // compute subcut
423 if( cut < DBL_MAX ) subcut = minSubRange*cut;
424 if(pcuts->GetProductionCut(idx) < maxCutInRange) {
425 G4double tcutmax =
426 theCoupleTable->ConvertRangeToEnergy(secondaryParticle,
427 material,maxSubCutInRange);
428 if(tcutmax < subcut) subcut = tcutmax;
429 }
430 }
431
432 G4int nm = setOfRegionModels[reg]->NumberOfModels();
433 for(G4int j=0; j<nm; j++) {
434
435 G4VEmModel* model = models[setOfRegionModels[reg]->ModelIndex(j)];
436
437 G4double tcutmin = model->MinEnergyCut(particle, couple);
438
439 if(cut < tcutmin) cut = tcutmin;
440 if(subcut < tcutmin) subcut = tcutmin;
441 if(1 < verboseLevel) {
442 G4cout << "The model # " << j
443 << "; tcutmin(MeV)= " << tcutmin/MeV
444 << "; tcut(MeV)= " << cut/MeV
445 << "; tsubcut(MeV)= " << subcut/MeV
446 << " for " << particle->GetParticleName()
447 << G4endl;
448 }
449 }
450 theCuts.push_back(cut);
451 theSubCuts.push_back(subcut);
452 }
453
454 for(G4int jj=0; jj<nEmModels; jj++) {
455 models[jj]->Initialise(particle, theCuts);
456 if(flucModels[jj]) flucModels[jj]->InitialiseMe(particle);
457 }
458
459
460 if(1 < verboseLevel) {
461 G4cout << "G4EmModelManager for " << particle->GetParticleName()
462 << " is initialised; nRegions= " << nRegions
463 << G4endl;
464 }
465
466 return &theCuts;
467}
468
469//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
470
471void G4EmModelManager::FillDEDXVector(G4PhysicsVector* aVector,
472 const G4MaterialCutsCouple* couple,
473 G4EmTableType tType)
474{
475
476 // vectors to provide continues dE/dx
477 G4DataVector factor;
478 G4DataVector dedxLow;
479 G4DataVector dedxHigh;
480
481 G4double e;
482
483 size_t i = couple->GetIndex();
484 G4double cut = theCuts[i];
485 G4double subcut = 0.0;
486
487 if(fTotal == tType) cut = DBL_MAX;
488 else if(fSubRestricted == tType) subcut = theSubCuts[i];
489
490 if(1 < verboseLevel) {
491 G4cout << "G4EmModelManager::FillDEDXVector() for "
492 << couple->GetMaterial()->GetName()
493 << " Ecut(MeV)= " << cut
494 << " Esubcut(MeV)= " << subcut
495 << " Type " << tType
496 << " for " << particle->GetParticleName()
497 << G4endl;
498 }
499
500 G4int reg = idxOfRegionModels[i];
501 const G4RegionModels* regModels = setOfRegionModels[reg];
502 G4int nmod = regModels->NumberOfModels();
503 factor.resize(nmod);
504 dedxLow.resize(nmod);
505 dedxHigh.resize(nmod);
506
507 if(1 < verboseLevel) {
508 G4cout << "There are " << nmod << " models for "
509 << couple->GetMaterial()->GetName()
510 << " at the region #" << reg
511 << G4endl;
512 }
513
514
515 // calculate factors to provide continuity of energy loss
516 factor[0] = 1.0;
517 G4int j;
518
519 G4int totBinsLoss = aVector->GetVectorLength();
520
521 dedxLow[0] = 0.0;
522
523 e = upperEkin[regModels->ModelIndex(0)];
524 G4VEmModel* model = models[regModels->ModelIndex(0)];
525 dedxHigh[0] = 0.0;
526 if(model && cut > subcut) {
527 dedxHigh[0] = model->ComputeDEDX(couple,particle,e,cut);
528 if(subcut > 0.0)
529 dedxHigh[0] -= model->ComputeDEDX(couple,particle,e,subcut);
530 }
531 if(nmod > 1) {
532 for(j=1; j<nmod; j++) {
533
534 e = upperEkin[regModels->ModelIndex(j-1)];
535 G4int idx = regModels->ModelIndex(j);
536
537 dedxLow[j] = models[idx]->ComputeDEDX(couple,particle,e,cut);
538 if(subcut > 0.0)
539 dedxLow[j] -= models[idx]->ComputeDEDX(couple,particle,e,subcut);
540 if(subcut == cut) dedxLow[j] = 0.0;
541
542 e = upperEkin[idx];
543 dedxHigh[j] = models[idx]->ComputeDEDX(couple,particle,e,cut);
544 if(subcut > 0.0)
545 dedxHigh[j] -= models[idx]->ComputeDEDX(couple,particle,e,subcut);
546 if(subcut == cut) dedxHigh[j] = 0.0;
547 }
548
549 for(j=1; j<nmod; j++) {
550 if(dedxLow[j] > 0.0) factor[j] = (dedxHigh[j-1]/dedxLow[j] - 1.0);
551 else factor[j] = 0.0;
552 }
553
554 if(2 < verboseLevel) {
555 G4cout << "Loop over " << totBinsLoss << " bins start " << G4endl;
556 }
557 }
558
559 // Calculate energy losses vector
560 for(j=0; j<totBinsLoss; j++) {
561
562 G4double e = aVector->GetLowEdgeEnergy(j);
563 G4double fac = 1.0;
564
565 // Choose a model of energy losses
566 G4int k = 0;
567 if (nmod > 1 && e > upperEkin[regModels->ModelIndex(0)]) {
568 do {
569 k++;
570 fac *= (1.0 + factor[k]*upperEkin[regModels->ModelIndex(k-1)]/e);
571 } while (k<nmod-1 && e > upperEkin[regModels->ModelIndex(k)] );
572 }
573
574 model = models[regModels->ModelIndex(k)];
575 G4double dedx = 0.0;
576 G4double dedx0 = 0.0;
577 if(model && cut > subcut) {
578 dedx = model->ComputeDEDX(couple,particle,e,cut);
579 dedx0 = dedx;
580 if(subcut > 0.0) dedx -= model->ComputeDEDX(couple,particle,e,subcut);
581 dedx *= fac;
582 }
583
584 if(dedx < 0.0) dedx = 0.0;
585 if(2 < verboseLevel) {
586 G4cout << "Material= " << couple->GetMaterial()->GetName()
587 << " E(MeV)= " << e/MeV
588 << " dEdx(MeV/mm)= " << dedx*mm/MeV
589 << " dEdx0(MeV/mm)= " << dedx0*mm/MeV
590 << " fac= " << fac
591 << G4endl;
592 }
593 aVector->PutValue(j, dedx);
594 }
595}
596
597//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
598
599void G4EmModelManager::FillLambdaVector(G4PhysicsVector* aVector,
600 const G4MaterialCutsCouple* couple,
601 G4bool startFromNull,
602 G4EmTableType tType)
603{
604 // vectors to provide continues cross section
605 G4DataVector factor;
606 G4DataVector sigmaLow;
607 G4DataVector sigmaHigh;
608
609 G4double e;
610
611 size_t i = couple->GetIndex();
612 G4double cut = theCuts[i];
613 G4double tmax = DBL_MAX;
614 if (fSubRestricted == tType) {
615 tmax = cut;
616 cut = theSubCuts[i];
617 }
618
619 if(1 < verboseLevel) {
620 G4cout << "G4EmModelManager::FillLambdaVector() for particle "
621 << particle->GetParticleName()
622 << " in " << couple->GetMaterial()->GetName()
623 << " Ecut(MeV)= " << cut
624 << " Emax(MeV)= " << tmax
625 << " Type " << tType
626 << G4endl;
627 }
628
629 G4int reg = idxOfRegionModels[i];
630 const G4RegionModels* regModels = setOfRegionModels[reg];
631 G4int nmod = regModels->NumberOfModels();
632 factor.resize(nmod);
633 sigmaLow.resize(nmod);
634 sigmaHigh.resize(nmod);
635
636 if(2 < verboseLevel) {
637 G4cout << "There are " << nmod << " models for "
638 << couple->GetMaterial()->GetName() << G4endl;
639 }
640
641 // calculate factors to provide continuity of energy loss
642 factor[0] = 1.0;
643 G4int j;
644 G4int totBinsLambda = aVector->GetVectorLength();
645
646 sigmaLow[0] = 0.0;
647
648 e = upperEkin[regModels->ModelIndex(0)];
649 G4VEmModel* model = models[regModels->ModelIndex(0)];
650 sigmaHigh[0] = 0.0;
651 if(model) sigmaHigh[0] = model->CrossSection(couple,particle,e,cut,tmax);
652
653 if(2 < verboseLevel) {
654 G4cout << "### For material " << couple->GetMaterial()->GetName()
655 << " " << nmod
656 << " models"
657 << " Ecut(MeV)= " << cut/MeV
658 << " Emax(MeV)= " << e/MeV
659 << " nbins= " << totBinsLambda
660 << G4endl;
661 G4cout << " model #0 eUp= " << e
662 << " sigmaUp= " << sigmaHigh[0] << G4endl;
663 }
664
665
666 if(nmod > 1) {
667
668 for(j=1; j<nmod; j++) {
669
670 e = upperEkin[regModels->ModelIndex(j-1)];
671 sigmaLow[j] =
672 models[regModels->ModelIndex(j)]->CrossSection(couple,particle,e,cut,tmax);
673 e = upperEkin[regModels->ModelIndex(j)];
674 sigmaHigh[j] =
675 models[regModels->ModelIndex(j)]->CrossSection(couple,particle,e,cut,tmax);
676 if(1 < verboseLevel) {
677 G4cout << " model #" << j << " eUp= " << e
678 << " sigmaUp= " << sigmaHigh[j]
679 << " sigmaDown= " << sigmaLow[j]
680 << G4endl;
681 }
682 }
683 for(j=1; j<nmod; j++) {
684 if(sigmaLow[j] > 0.0) factor[j] = (sigmaHigh[j-1]/sigmaLow[j] - 1.0);
685 else factor[j] = 0.0;
686 }
687 }
688
689 // Calculate lambda vector
690 for(j=0; j<totBinsLambda; j++) {
691
692 e = aVector->GetLowEdgeEnergy(j);
693
694 // Choose a model of energy losses
695 G4int k = 0;
696 G4double fac = 1.0;
697 if (nmod > 1 && e > upperEkin[regModels->ModelIndex(0)]) {
698 do {
699 k++;
700 fac *= (1.0 + factor[k]*upperEkin[regModels->ModelIndex(k-1)]/e);
701 } while (k<nmod-1 && e > upperEkin[regModels->ModelIndex(k)] );
702 }
703
704 model = models[regModels->ModelIndex(k)];
705 G4double cross = 0.0;
706 if(model) cross = model->CrossSection(couple,particle,e,cut,tmax)*fac;
707 if(j==0 && startFromNull) cross = 0.0;
708
709 if(2 < verboseLevel) {
710 G4cout << "FillLambdaVector: " << j << ". e(MeV)= " << e/MeV
711 << " cross(1/mm)= " << cross*mm
712 << " fac= " << fac << " k= " << k
713 << " model= " << regModels->ModelIndex(k)
714 << G4endl;
715 }
716 if(cross < 0.0) cross = 0.0;
717
718 aVector->PutValue(j, cross);
719 }
720}
721
722//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
Note: See TracBrowser for help on using the repository browser.