source: trunk/source/processes/electromagnetic/utils/src/G4VEmProcess.cc@ 1181

Last change on this file since 1181 was 1055, checked in by garnier, 17 years ago

maj sur la beta de geant 4.9.3
File size: 21.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: G4VEmProcess.cc,v 1.66 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: G4VEmProcess
35//
36// Author: Vladimir Ivanchenko on base of Laszlo Urban code
37//
38// Creation date: 01.10.2003
39//
40// Modifications:
41// 30-06-04 make it to be pure discrete process (V.Ivanchenko)
42// 30-09-08 optimise integral option (V.Ivanchenko)
43// 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivanchenko)
44// 11-03-05 Shift verbose level by 1, add applyCuts and killPrimary flags (VI)
45// 14-03-05 Update logic PostStepDoIt (V.Ivanchenko)
46// 08-04-05 Major optimisation of internal interfaces (V.Ivanchenko)
47// 18-04-05 Use G4ParticleChangeForGamma (V.Ivanchenko)
48// 25-07-05 Add protection: integral mode only for charged particles (VI)
49// 04-09-05 default lambdaFactor 0.8 (V.Ivanchenko)
50// 11-01-06 add A to parameters of ComputeCrossSectionPerAtom (VI)
51// 12-09-06 add SetModel() (mma)
52// 12-04-07 remove double call to Clear model manager (V.Ivanchenko)
53// 27-10-07 Virtual functions moved to source (V.Ivanchenko)
54//
55// Class Description:
56//
57
58// -------------------------------------------------------------------
59//
60//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
61//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
62
63#include "G4VEmProcess.hh"
64#include "G4LossTableManager.hh"
65#include "G4Step.hh"
66#include "G4ParticleDefinition.hh"
67#include "G4VEmModel.hh"
68#include "G4DataVector.hh"
69#include "G4PhysicsTable.hh"
70#include "G4PhysicsVector.hh"
71#include "G4PhysicsLogVector.hh"
72#include "G4VParticleChange.hh"
73#include "G4ProductionCutsTable.hh"
74#include "G4Region.hh"
75#include "G4RegionStore.hh"
76#include "G4Gamma.hh"
77#include "G4Electron.hh"
78#include "G4Positron.hh"
79#include "G4PhysicsTableHelper.hh"
80#include "G4EmConfigurator.hh"
81
82//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
83
84G4VEmProcess::G4VEmProcess(const G4String& name, G4ProcessType type):
85 G4VDiscreteProcess(name, type),
86 secondaryParticle(0),
87 buildLambdaTable(true),
88 theLambdaTable(0),
89 theEnergyOfCrossSectionMax(0),
90 theCrossSectionMax(0),
91 integral(false),
92 applyCuts(false),
93 startFromNull(true),
94 useDeexcitation(false),
95 nDERegions(0),
96 idxDERegions(0),
97 currentModel(0),
98 particle(0),
99 currentCouple(0)
100{
101 SetVerboseLevel(1);
102
103 // Size of tables assuming spline
104 minKinEnergy = 0.1*keV;
105 maxKinEnergy = 10.0*TeV;
106 nLambdaBins = 77;
107
108 // default lambda factor
109 lambdaFactor = 0.8;
110
111 // default limit on polar angle
112 polarAngleLimit = 0.0;
113
114 // particle types
115 theGamma = G4Gamma::Gamma();
116 theElectron = G4Electron::Electron();
117 thePositron = G4Positron::Positron();
118
119 pParticleChange = &fParticleChange;
120 secParticles.reserve(5);
121
122 modelManager = new G4EmModelManager();
123 (G4LossTableManager::Instance())->Register(this);
124}
125
126//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
127
128G4VEmProcess::~G4VEmProcess()
129{
130 if(1 < verboseLevel)
131 G4cout << "G4VEmProcess destruct " << GetProcessName()
132 << G4endl;
133 Clear();
134 if(theLambdaTable) {
135 theLambdaTable->clearAndDestroy();
136 delete theLambdaTable;
137 }
138 delete modelManager;
139 (G4LossTableManager::Instance())->DeRegister(this);
140}
141
142//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
143
144void G4VEmProcess::Clear()
145{
146 delete [] theEnergyOfCrossSectionMax;
147 delete [] theCrossSectionMax;
148 delete [] idxDERegions;
149 theEnergyOfCrossSectionMax = 0;
150 theCrossSectionMax = 0;
151 idxDERegions = 0;
152 currentCouple = 0;
153 preStepLambda = 0.0;
154 mfpKinEnergy = DBL_MAX;
155 deRegions.clear();
156 nDERegions = 0;
157}
158
159//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
160
161void G4VEmProcess::AddEmModel(G4int order, G4VEmModel* p,
162 const G4Region* region)
163{
164 G4VEmFluctuationModel* fm = 0;
165 modelManager->AddEmModel(order, p, fm, region);
166 if(p) p->SetParticleChange(pParticleChange);
167}
168
169//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
170
171void G4VEmProcess::SetModel(G4VEmModel* p, G4int index)
172{
173 G4int n = emModels.size();
174 if(index >= n) for(G4int i=n; i<index+1; i++) {emModels.push_back(0);}
175 emModels[index] = p;
176}
177
178//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
179
180G4VEmModel* G4VEmProcess::Model(G4int index)
181{
182 G4VEmModel* p = 0;
183 if(index >= 0 && index < G4int(emModels.size())) p = emModels[index];
184 return p;
185}
186
187//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
188
189void G4VEmProcess::UpdateEmModel(const G4String& nam,
190 G4double emin, G4double emax)
191{
192 modelManager->UpdateEmModel(nam, emin, emax);
193}
194
195//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
196
197G4VEmModel* G4VEmProcess::GetModelByIndex(G4int idx, G4bool ver)
198{
199 return modelManager->GetModel(idx, ver);
200}
201
202//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
203
204void G4VEmProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
205{
206 if(!particle) particle = &part;
207 if(1 < verboseLevel) {
208 G4cout << "G4VEmProcess::PreparePhysicsTable() for "
209 << GetProcessName()
210 << " and particle " << part.GetParticleName()
211 << " local particle " << particle->GetParticleName()
212 << G4endl;
213 }
214
215 (G4LossTableManager::Instance())->EmConfigurator()->AddModels();
216
217 if(particle == &part) {
218 Clear();
219 InitialiseProcess(particle);
220 theCuts = modelManager->Initialise(particle,secondaryParticle,2.,verboseLevel);
221 const G4ProductionCutsTable* theCoupleTable=
222 G4ProductionCutsTable::GetProductionCutsTable();
223 theCutsGamma = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut);
224 theCutsElectron = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut);
225 theCutsPositron = theCoupleTable->GetEnergyCutsVector(idxG4PositronCut);
226 if(buildLambdaTable){
227 theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable);
228 }
229 }
230 // Sub Cutoff and Deexcitation
231 if (nDERegions>0) {
232
233 const G4ProductionCutsTable* theCoupleTable=
234 G4ProductionCutsTable::GetProductionCutsTable();
235 size_t numOfCouples = theCoupleTable->GetTableSize();
236
237 idxDERegions = new G4bool[numOfCouples];
238
239 for (size_t j=0; j<numOfCouples; j++) {
240
241 const G4MaterialCutsCouple* couple =
242 theCoupleTable->GetMaterialCutsCouple(j);
243 const G4ProductionCuts* pcuts = couple->GetProductionCuts();
244 G4bool reg = false;
245 for(G4int i=0; i<nDERegions; i++) {
246 if(deRegions[i]) {
247 if(pcuts == deRegions[i]->GetProductionCuts()) reg = true;
248 }
249 }
250 idxDERegions[j] = reg;
251 }
252 }
253 if (1 < verboseLevel && nDERegions>0) {
254 G4cout << " Deexcitation is activated for regions: " << G4endl;
255 for (G4int i=0; i<nDERegions; i++) {
256 const G4Region* r = deRegions[i];
257 G4cout << " " << r->GetName() << G4endl;
258 }
259 }
260}
261
262//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
263
264void G4VEmProcess::BuildPhysicsTable(const G4ParticleDefinition& part)
265{
266 if(1 < verboseLevel) {
267 G4cout << "G4VEmProcess::BuildPhysicsTable() for "
268 << GetProcessName()
269 << " and particle " << part.GetParticleName()
270 << " buildLambdaTable= " << buildLambdaTable
271 << G4endl;
272 }
273
274 if(buildLambdaTable) {
275 BuildLambdaTable();
276 FindLambdaMax();
277 }
278 if(0 < verboseLevel) PrintInfoDefinition();
279
280 if(1 < verboseLevel) {
281 G4cout << "G4VEmProcess::BuildPhysicsTable() done for "
282 << GetProcessName()
283 << " and particle " << part.GetParticleName()
284 << G4endl;
285 }
286}
287
288//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
289
290void G4VEmProcess::BuildLambdaTable()
291{
292 if(1 < verboseLevel) {
293 G4cout << "G4EmProcess::BuildLambdaTable() for process "
294 << GetProcessName() << " and particle "
295 << particle->GetParticleName()
296 << G4endl;
297 }
298
299 // Access to materials
300 const G4ProductionCutsTable* theCoupleTable=
301 G4ProductionCutsTable::GetProductionCutsTable();
302 size_t numOfCouples = theCoupleTable->GetTableSize();
303 for(size_t i=0; i<numOfCouples; i++) {
304
305 if (theLambdaTable->GetFlag(i)) {
306
307 // create physics vector and fill it
308 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
309 G4PhysicsVector* aVector = LambdaPhysicsVector(couple);
310 modelManager->FillLambdaVector(aVector, couple, startFromNull);
311 G4PhysicsTableHelper::SetPhysicsVector(theLambdaTable, i, aVector);
312 }
313 }
314
315 if(1 < verboseLevel) {
316 G4cout << "Lambda table is built for "
317 << particle->GetParticleName()
318 << G4endl;
319 if(2 < verboseLevel) {
320 G4cout << *theLambdaTable << G4endl;
321 }
322 }
323}
324
325//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
326
327void G4VEmProcess::PrintInfoDefinition()
328{
329 if(verboseLevel > 0) {
330 G4cout << G4endl << GetProcessName() << ": for "
331 << particle->GetParticleName();
332 if(integral) G4cout << ", integral: 1 ";
333 if(applyCuts) G4cout << ", applyCuts: 1 ";
334 G4cout << " SubType= " << GetProcessSubType() << G4endl;
335 if(buildLambdaTable) {
336 G4cout << " Lambda tables from "
337 << G4BestUnit(minKinEnergy,"Energy")
338 << " to "
339 << G4BestUnit(maxKinEnergy,"Energy")
340 << " in " << nLambdaBins << " bins, spline: "
341 << (G4LossTableManager::Instance())->SplineFlag()
342 << G4endl;
343 }
344 PrintInfo();
345 modelManager->DumpModelList(verboseLevel);
346 }
347
348 if(verboseLevel > 2 && buildLambdaTable) {
349 G4cout << " LambdaTable address= " << theLambdaTable << G4endl;
350 if(theLambdaTable) G4cout << (*theLambdaTable) << G4endl;
351 }
352}
353
354//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
355
356G4double G4VEmProcess::PostStepGetPhysicalInteractionLength(
357 const G4Track& track,
358 G4double previousStepSize,
359 G4ForceCondition* condition)
360{
361 // condition is set to "Not Forced"
362 *condition = NotForced;
363 G4double x = DBL_MAX;
364 if(previousStepSize <= DBL_MIN) theNumberOfInteractionLengthLeft = -1.0;
365 InitialiseStep(track);
366
367 if(preStepKinEnergy < mfpKinEnergy) {
368 if (integral) ComputeIntegralLambda(preStepKinEnergy);
369 else preStepLambda = GetCurrentLambda(preStepKinEnergy);
370 if(preStepLambda <= DBL_MIN) mfpKinEnergy = 0.0;
371 }
372
373 // non-zero cross section
374 if(preStepLambda > DBL_MIN) {
375 if (theNumberOfInteractionLengthLeft < 0.0) {
376 // beggining of tracking (or just after DoIt of this process)
377 ResetNumberOfInteractionLengthLeft();
378 } else if(currentInteractionLength < DBL_MAX) {
379 // subtract NumberOfInteractionLengthLeft
380 SubtractNumberOfInteractionLengthLeft(previousStepSize);
381 if(theNumberOfInteractionLengthLeft < 0.)
382 theNumberOfInteractionLengthLeft = perMillion;
383 }
384
385 // get mean free path and step limit
386 currentInteractionLength = 1.0/preStepLambda;
387 x = theNumberOfInteractionLengthLeft * currentInteractionLength;
388#ifdef G4VERBOSE
389 if (verboseLevel>2){
390 G4cout << "G4VEmProcess::PostStepGetPhysicalInteractionLength ";
391 G4cout << "[ " << GetProcessName() << "]" << G4endl;
392 G4cout << " for " << particle->GetParticleName()
393 << " in Material " << currentMaterial->GetName()
394 << " Ekin(MeV)= " << preStepKinEnergy/MeV
395 <<G4endl;
396 G4cout << "MeanFreePath = " << currentInteractionLength/cm << "[cm]"
397 << "InteractionLength= " << x/cm <<"[cm] " <<G4endl;
398 }
399#endif
400
401 // zero cross section case
402 } else {
403 if(theNumberOfInteractionLengthLeft > DBL_MIN &&
404 currentInteractionLength < DBL_MAX) {
405
406 // subtract NumberOfInteractionLengthLeft
407 SubtractNumberOfInteractionLengthLeft(previousStepSize);
408 if(theNumberOfInteractionLengthLeft < 0.)
409 theNumberOfInteractionLengthLeft = perMillion;
410 }
411 currentInteractionLength = DBL_MAX;
412 }
413 return x;
414}
415
416//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
417
418G4VParticleChange* G4VEmProcess::PostStepDoIt(const G4Track& track,
419 const G4Step&)
420{
421 fParticleChange.InitializeForPostStep(track);
422
423 // Do not make anything if particle is stopped, the annihilation then
424 // should be performed by the AtRestDoIt!
425 if (track.GetTrackStatus() == fStopButAlive) return &fParticleChange;
426
427 G4double finalT = track.GetKineticEnergy();
428
429 // Integral approach
430 if (integral) {
431 G4double lx = GetLambda(finalT, currentCouple);
432 if(preStepLambda<lx && 1 < verboseLevel) {
433 G4cout << "WARING: for " << particle->GetParticleName()
434 << " and " << GetProcessName()
435 << " E(MeV)= " << finalT/MeV
436 << " preLambda= " << preStepLambda << " < " << lx << " (postLambda) "
437 << G4endl;
438 }
439
440 if(preStepLambda*G4UniformRand() > lx) {
441 ClearNumberOfInteractionLengthLeft();
442 return &fParticleChange;
443 }
444 }
445
446 SelectModel(finalT);
447 if(useDeexcitation) {
448 currentModel->SetDeexcitationFlag(idxDERegions[currentMaterialIndex]);
449 }
450 /*
451 if(0 < verboseLevel) {
452 G4cout << "G4VEmProcess::PostStepDoIt: Sample secondary; E= "
453 << finalT/MeV
454 << " MeV; model= (" << currentModel->LowEnergyLimit()
455 << ", " << currentModel->HighEnergyLimit() << ")"
456 << G4endl;
457 }
458 */
459
460
461 // sample secondaries
462 secParticles.clear();
463 currentModel->SampleSecondaries(&secParticles,
464 currentCouple,
465 track.GetDynamicParticle(),
466 (*theCuts)[currentMaterialIndex]);
467
468 // save secondaries
469 G4int num = secParticles.size();
470 if(num > 0) {
471
472 fParticleChange.SetNumberOfSecondaries(num);
473 G4double edep = fParticleChange.GetLocalEnergyDeposit();
474
475 for (G4int i=0; i<num; i++) {
476 G4DynamicParticle* dp = secParticles[i];
477 const G4ParticleDefinition* p = dp->GetDefinition();
478 G4double e = dp->GetKineticEnergy();
479 G4bool good = true;
480 if(applyCuts) {
481 if (p == theGamma) {
482 if (e < (*theCutsGamma)[currentMaterialIndex]) good = false;
483
484 } else if (p == theElectron) {
485 if (e < (*theCutsElectron)[currentMaterialIndex]) good = false;
486
487 } else if (p == thePositron) {
488 if (electron_mass_c2 < (*theCutsGamma)[currentMaterialIndex] &&
489 e < (*theCutsPositron)[currentMaterialIndex]) {
490 good = false;
491 e += 2.0*electron_mass_c2;
492 }
493 }
494 if(!good) {
495 delete dp;
496 edep += e;
497 }
498 }
499 if (good) fParticleChange.AddSecondary(dp);
500 }
501 fParticleChange.ProposeLocalEnergyDeposit(edep);
502 }
503
504 ClearNumberOfInteractionLengthLeft();
505 return &fParticleChange;
506}
507
508//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
509
510G4bool G4VEmProcess::StorePhysicsTable(const G4ParticleDefinition* part,
511 const G4String& directory,
512 G4bool ascii)
513{
514 G4bool yes = true;
515
516 if ( theLambdaTable && part == particle) {
517 const G4String name =
518 GetPhysicsTableFileName(part,directory,"Lambda",ascii);
519 yes = theLambdaTable->StorePhysicsTable(name,ascii);
520
521 if ( yes ) {
522 G4cout << "Physics tables are stored for " << particle->GetParticleName()
523 << " and process " << GetProcessName()
524 << " in the directory <" << directory
525 << "> " << G4endl;
526 } else {
527 G4cout << "Fail to store Physics Tables for "
528 << particle->GetParticleName()
529 << " and process " << GetProcessName()
530 << " in the directory <" << directory
531 << "> " << G4endl;
532 }
533 }
534 return yes;
535}
536
537//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
538
539G4bool G4VEmProcess::RetrievePhysicsTable(const G4ParticleDefinition* part,
540 const G4String& directory,
541 G4bool ascii)
542{
543 if(1 < verboseLevel) {
544 G4cout << "G4VEmProcess::RetrievePhysicsTable() for "
545 << part->GetParticleName() << " and process "
546 << GetProcessName() << G4endl;
547 }
548 G4bool yes = true;
549
550 if(!buildLambdaTable || particle != part) return yes;
551
552 const G4String particleName = part->GetParticleName();
553 G4String filename;
554
555 filename = GetPhysicsTableFileName(part,directory,"Lambda",ascii);
556 yes = G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTable,
557 filename,ascii);
558 if ( yes ) {
559 if (0 < verboseLevel) {
560 G4cout << "Lambda table for " << particleName << " is Retrieved from <"
561 << filename << ">"
562 << G4endl;
563 }
564 if((G4LossTableManager::Instance())->SplineFlag()) {
565 size_t n = theLambdaTable->length();
566 for(size_t i=0; i<n; i++) {(* theLambdaTable)[i]->SetSpline(true);}
567 }
568 } else {
569 if (1 < verboseLevel) {
570 G4cout << "Lambda table for " << particleName << " in file <"
571 << filename << "> is not exist"
572 << G4endl;
573 }
574 }
575
576 return yes;
577}
578
579//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
580
581void G4VEmProcess::ActivateDeexcitation(G4bool val, const G4Region* r)
582{
583 G4RegionStore* regionStore = G4RegionStore::GetInstance();
584 const G4Region* reg = r;
585 if (!reg) {reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);}
586
587 // the region is in the list
588 if (nDERegions) {
589 for (G4int i=0; i<nDERegions; i++) {
590 if (reg == deRegions[i]) {
591 if(!val) deRegions[i] = 0;
592 return;
593 }
594 }
595 }
596
597 // new region
598 if(val) {
599 useDeexcitation = true;
600 deRegions.push_back(reg);
601 nDERegions++;
602 } else {
603 useDeexcitation = false;
604 }
605}
606
607//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
608
609G4double G4VEmProcess::CrossSectionPerVolume(G4double kineticEnergy,
610 const G4MaterialCutsCouple* couple)
611{
612 // Cross section per atom is calculated
613 DefineMaterial(couple);
614 G4double cross = 0.0;
615 G4bool b;
616 if(theLambdaTable) {
617 cross = (((*theLambdaTable)[currentMaterialIndex])->
618 GetValue(kineticEnergy, b));
619 } else {
620 SelectModel(kineticEnergy);
621 cross = currentModel->CrossSectionPerVolume(currentMaterial,
622 particle,kineticEnergy);
623 }
624
625 return cross;
626}
627
628//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
629
630G4double G4VEmProcess::GetMeanFreePath(const G4Track& track,
631 G4double,
632 G4ForceCondition* condition)
633{
634 *condition = NotForced;
635 return G4VEmProcess::MeanFreePath(track);
636}
637
638//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
639
640void G4VEmProcess::FindLambdaMax()
641{
642 if(1 < verboseLevel) {
643 G4cout << "### G4VEmProcess::FindLambdaMax: "
644 << particle->GetParticleName()
645 << " and process " << GetProcessName() << G4endl;
646 }
647 size_t n = theLambdaTable->length();
648 G4PhysicsVector* pv = (*theLambdaTable)[0];
649 G4double e, s, emax, smax;
650 theEnergyOfCrossSectionMax = new G4double [n];
651 theCrossSectionMax = new G4double [n];
652 G4bool b;
653
654 for (size_t i=0; i<n; i++) {
655 pv = (*theLambdaTable)[i];
656 emax = DBL_MAX;
657 smax = 0.0;
658 if(pv) {
659 size_t nb = pv->GetVectorLength();
660 emax = pv->GetLowEdgeEnergy(nb);
661 smax = 0.0;
662 for (size_t j=0; j<nb; j++) {
663 e = pv->GetLowEdgeEnergy(j);
664 s = pv->GetValue(e,b);
665 if(s > smax) {
666 smax = s;
667 emax = e;
668 }
669 }
670 }
671 theEnergyOfCrossSectionMax[i] = emax;
672 theCrossSectionMax[i] = smax;
673 if(2 < verboseLevel) {
674 G4cout << "For " << particle->GetParticleName()
675 << " Max CS at i= " << i << " emax(MeV)= " << emax/MeV
676 << " lambda= " << smax << G4endl;
677 }
678 }
679}
680
681//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
682
683G4PhysicsVector* G4VEmProcess::LambdaPhysicsVector(const G4MaterialCutsCouple*)
684{
685 G4PhysicsVector* v =
686 new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nLambdaBins);
687 v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
688 return v;
689}
690
691//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
Note: See TracBrowser for help on using the repository browser.