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

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

update to geant4.9.2
File size: 18.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.60 2008/10/17 14:46:16 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-02 $
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
81//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
82
83G4VEmProcess::G4VEmProcess(const G4String& name, G4ProcessType type):
84 G4VDiscreteProcess(name, type),
85 secondaryParticle(0),
86 buildLambdaTable(true),
87 theLambdaTable(0),
88 theEnergyOfCrossSectionMax(0),
89 theCrossSectionMax(0),
90 integral(false),
91 applyCuts(false),
92 startFromNull(true),
93 nRegions(0),
94 selectedModel(0),
95 particle(0),
96 currentCouple(0)
97{
98 SetVerboseLevel(1);
99
100 // Size of tables assuming spline
101 minKinEnergy = 0.1*keV;
102 maxKinEnergy = 100.0*TeV;
103 nLambdaBins = 84;
104
105 // default lambda factor
106 lambdaFactor = 0.8;
107
108 // default limit on polar angle
109 polarAngleLimit = 0.0;
110
111 // particle types
112 theGamma = G4Gamma::Gamma();
113 theElectron = G4Electron::Electron();
114 thePositron = G4Positron::Positron();
115
116 pParticleChange = &fParticleChange;
117 secParticles.reserve(5);
118
119 modelManager = new G4EmModelManager();
120 (G4LossTableManager::Instance())->Register(this);
121}
122
123//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
124
125G4VEmProcess::~G4VEmProcess()
126{
127 if(1 < verboseLevel)
128 G4cout << "G4VEmProcess destruct " << GetProcessName()
129 << G4endl;
130 Clear();
131 if(theLambdaTable) {
132 theLambdaTable->clearAndDestroy();
133 delete theLambdaTable;
134 }
135 delete modelManager;
136 (G4LossTableManager::Instance())->DeRegister(this);
137}
138
139//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
140
141void G4VEmProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
142{
143 if(!particle) particle = &part;
144 if(1 < verboseLevel) {
145 G4cout << "G4VEmProcess::PreparePhysicsTable() for "
146 << GetProcessName()
147 << " and particle " << part.GetParticleName()
148 << " local particle " << particle->GetParticleName()
149 << G4endl;
150 }
151
152 if(particle == &part) {
153 Clear();
154 InitialiseProcess(particle);
155 theCuts = modelManager->Initialise(particle,secondaryParticle,2.,verboseLevel);
156 const G4ProductionCutsTable* theCoupleTable=
157 G4ProductionCutsTable::GetProductionCutsTable();
158 theCutsGamma = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut);
159 theCutsElectron = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut);
160 theCutsPositron = theCoupleTable->GetEnergyCutsVector(idxG4PositronCut);
161 if(buildLambdaTable)
162 theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable);
163 }
164}
165
166//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
167
168void G4VEmProcess::Clear()
169{
170 if(theEnergyOfCrossSectionMax) delete [] theEnergyOfCrossSectionMax;
171 if(theCrossSectionMax) delete [] theCrossSectionMax;
172 theEnergyOfCrossSectionMax = 0;
173 theCrossSectionMax = 0;
174 currentCouple = 0;
175 preStepLambda = 0.0;
176 mfpKinEnergy = DBL_MAX;
177}
178
179//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
180
181void G4VEmProcess::BuildPhysicsTable(const G4ParticleDefinition& part)
182{
183 if(1 < verboseLevel) {
184 G4cout << "G4VEmProcess::BuildPhysicsTable() for "
185 << GetProcessName()
186 << " and particle " << part.GetParticleName()
187 << " buildLambdaTable= " << buildLambdaTable
188 << G4endl;
189 }
190
191 if(buildLambdaTable) {
192 BuildLambdaTable();
193 FindLambdaMax();
194 }
195 if(0 < verboseLevel) PrintInfoDefinition();
196
197 if(1 < verboseLevel) {
198 G4cout << "G4VEmProcess::BuildPhysicsTable() done for "
199 << GetProcessName()
200 << " and particle " << part.GetParticleName()
201 << G4endl;
202 }
203}
204
205//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
206
207void G4VEmProcess::BuildLambdaTable()
208{
209 if(1 < verboseLevel) {
210 G4cout << "G4EmProcess::BuildLambdaTable() for process "
211 << GetProcessName() << " and particle "
212 << particle->GetParticleName()
213 << G4endl;
214 }
215
216 // Access to materials
217 const G4ProductionCutsTable* theCoupleTable=
218 G4ProductionCutsTable::GetProductionCutsTable();
219 size_t numOfCouples = theCoupleTable->GetTableSize();
220 for(size_t i=0; i<numOfCouples; i++) {
221
222 if (theLambdaTable->GetFlag(i)) {
223
224 // create physics vector and fill it
225 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
226 G4PhysicsVector* aVector = LambdaPhysicsVector(couple);
227 modelManager->FillLambdaVector(aVector, couple, startFromNull);
228 G4PhysicsTableHelper::SetPhysicsVector(theLambdaTable, i, aVector);
229 }
230 }
231
232 if(1 < verboseLevel) {
233 G4cout << "Lambda table is built for "
234 << particle->GetParticleName()
235 << G4endl;
236 if(2 < verboseLevel) {
237 G4cout << *theLambdaTable << G4endl;
238 }
239 }
240}
241
242//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
243
244G4double G4VEmProcess::PostStepGetPhysicalInteractionLength(
245 const G4Track& track,
246 G4double previousStepSize,
247 G4ForceCondition* condition)
248{
249 // condition is set to "Not Forced"
250 *condition = NotForced;
251 G4double x = DBL_MAX;
252 if(previousStepSize <= DBL_MIN) theNumberOfInteractionLengthLeft = -1.0;
253 InitialiseStep(track);
254
255 if(preStepKinEnergy < mfpKinEnergy) {
256 if (integral) ComputeIntegralLambda(preStepKinEnergy);
257 else preStepLambda = GetCurrentLambda(preStepKinEnergy);
258 if(preStepLambda <= DBL_MIN) mfpKinEnergy = 0.0;
259 }
260
261 // non-zero cross section
262 if(preStepLambda > DBL_MIN) {
263 if (theNumberOfInteractionLengthLeft < 0.0) {
264 // beggining of tracking (or just after DoIt of this process)
265 ResetNumberOfInteractionLengthLeft();
266 } else if(currentInteractionLength < DBL_MAX) {
267 // subtract NumberOfInteractionLengthLeft
268 SubtractNumberOfInteractionLengthLeft(previousStepSize);
269 if(theNumberOfInteractionLengthLeft < 0.)
270 theNumberOfInteractionLengthLeft = perMillion;
271 }
272
273 // get mean free path and step limit
274 currentInteractionLength = 1.0/preStepLambda;
275 x = theNumberOfInteractionLengthLeft * currentInteractionLength;
276#ifdef G4VERBOSE
277 if (verboseLevel>2){
278 G4cout << "G4VEmProcess::PostStepGetPhysicalInteractionLength ";
279 G4cout << "[ " << GetProcessName() << "]" << G4endl;
280 G4cout << " for " << particle->GetParticleName()
281 << " in Material " << currentMaterial->GetName()
282 << " Ekin(MeV)= " << preStepKinEnergy/MeV
283 <<G4endl;
284 G4cout << "MeanFreePath = " << currentInteractionLength/cm << "[cm]"
285 << "InteractionLength= " << x/cm <<"[cm] " <<G4endl;
286 }
287#endif
288
289 // zero cross section case
290 } else {
291 if(theNumberOfInteractionLengthLeft > DBL_MIN &&
292 currentInteractionLength < DBL_MAX) {
293
294 // subtract NumberOfInteractionLengthLeft
295 SubtractNumberOfInteractionLengthLeft(previousStepSize);
296 if(theNumberOfInteractionLengthLeft < 0.)
297 theNumberOfInteractionLengthLeft = perMillion;
298 }
299 currentInteractionLength = DBL_MAX;
300 }
301 return x;
302}
303
304//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
305
306G4double G4VEmProcess::GetMeanFreePath(const G4Track& track,
307 G4double,
308 G4ForceCondition* condition)
309{
310 *condition = NotForced;
311 return G4VEmProcess::MeanFreePath(track);
312}
313
314//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
315
316G4VParticleChange* G4VEmProcess::PostStepDoIt(const G4Track& track,
317 const G4Step&)
318{
319 fParticleChange.InitializeForPostStep(track);
320
321 // Do not make anything if particle is stopped, the annihilation then
322 // should be performed by the AtRestDoIt!
323 if (track.GetTrackStatus() == fStopButAlive) return &fParticleChange;
324
325 G4double finalT = track.GetKineticEnergy();
326
327 // Integral approach
328 if (integral) {
329 G4double lx = GetLambda(finalT, currentCouple);
330 if(preStepLambda<lx && 1 < verboseLevel) {
331 G4cout << "WARING: for " << particle->GetParticleName()
332 << " and " << GetProcessName()
333 << " E(MeV)= " << finalT/MeV
334 << " preLambda= " << preStepLambda << " < " << lx << " (postLambda) "
335 << G4endl;
336 }
337
338 if(preStepLambda*G4UniformRand() > lx) {
339 ClearNumberOfInteractionLengthLeft();
340 return &fParticleChange;
341 }
342 }
343
344 G4VEmModel* currentModel = SelectModel(finalT);
345
346 /*
347 if(0 < verboseLevel) {
348 G4cout << "G4VEmProcess::PostStepDoIt: Sample secondary; E= "
349 << finalT/MeV
350 << " MeV; model= (" << currentModel->LowEnergyLimit()
351 << ", " << currentModel->HighEnergyLimit() << ")"
352 << G4endl;
353 }
354 */
355
356
357 // sample secondaries
358 secParticles.clear();
359 currentModel->SampleSecondaries(&secParticles,
360 currentCouple,
361 track.GetDynamicParticle(),
362 (*theCuts)[currentMaterialIndex]);
363
364 // save secondaries
365 G4int num = secParticles.size();
366 if(num > 0) {
367
368 fParticleChange.SetNumberOfSecondaries(num);
369 G4double edep = fParticleChange.GetLocalEnergyDeposit();
370
371 for (G4int i=0; i<num; i++) {
372 G4DynamicParticle* dp = secParticles[i];
373 const G4ParticleDefinition* p = dp->GetDefinition();
374 G4double e = dp->GetKineticEnergy();
375 G4bool good = true;
376 if(applyCuts) {
377 if (p == theGamma) {
378 if (e < (*theCutsGamma)[currentMaterialIndex]) good = false;
379
380 } else if (p == theElectron) {
381 if (e < (*theCutsElectron)[currentMaterialIndex]) good = false;
382
383 } else if (p == thePositron) {
384 if (electron_mass_c2 < (*theCutsGamma)[currentMaterialIndex] &&
385 e < (*theCutsPositron)[currentMaterialIndex]) {
386 good = false;
387 e += 2.0*electron_mass_c2;
388 }
389 }
390 if(!good) {
391 delete dp;
392 edep += e;
393 }
394 }
395 if (good) fParticleChange.AddSecondary(dp);
396 }
397 fParticleChange.ProposeLocalEnergyDeposit(edep);
398 }
399
400 ClearNumberOfInteractionLengthLeft();
401 return &fParticleChange;
402}
403
404//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
405
406void G4VEmProcess::PrintInfoDefinition()
407{
408 if(verboseLevel > 0) {
409 G4cout << G4endl << GetProcessName() << ": for "
410 << particle->GetParticleName();
411 if(integral) G4cout << ", integral: 1 ";
412 if(applyCuts) G4cout << ", applyCuts: 1 ";
413 G4cout << " SubType= " << GetProcessSubType() << G4endl;
414 if(buildLambdaTable) {
415 G4cout << " Lambda tables from "
416 << G4BestUnit(minKinEnergy,"Energy")
417 << " to "
418 << G4BestUnit(maxKinEnergy,"Energy")
419 << " in " << nLambdaBins << " bins, spline: "
420 << (G4LossTableManager::Instance())->SplineFlag()
421 << G4endl;
422 }
423 PrintInfo();
424 modelManager->DumpModelList(verboseLevel);
425 }
426
427 if(verboseLevel > 2 && buildLambdaTable) {
428 G4cout << " LambdaTable address= " << theLambdaTable << G4endl;
429 if(theLambdaTable) G4cout << (*theLambdaTable) << G4endl;
430 }
431}
432
433//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
434
435G4double G4VEmProcess::CrossSectionPerVolume(G4double kineticEnergy,
436 const G4MaterialCutsCouple* couple)
437{
438 // Cross section per atom is calculated
439 DefineMaterial(couple);
440 G4double cross = 0.0;
441 G4bool b;
442 if(theLambdaTable) {
443 cross = (((*theLambdaTable)[currentMaterialIndex])->
444 GetValue(kineticEnergy, b));
445 } else {
446 G4VEmModel* model = SelectModel(kineticEnergy);
447 cross =
448 model->CrossSectionPerVolume(currentMaterial,particle,kineticEnergy);
449 }
450
451 return cross;
452}
453
454//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
455
456G4bool G4VEmProcess::StorePhysicsTable(const G4ParticleDefinition* part,
457 const G4String& directory,
458 G4bool ascii)
459{
460 G4bool yes = true;
461
462 if ( theLambdaTable && part == particle) {
463 const G4String name =
464 GetPhysicsTableFileName(part,directory,"Lambda",ascii);
465 yes = theLambdaTable->StorePhysicsTable(name,ascii);
466
467 if ( yes ) {
468 G4cout << "Physics tables are stored for " << particle->GetParticleName()
469 << " and process " << GetProcessName()
470 << " in the directory <" << directory
471 << "> " << G4endl;
472 } else {
473 G4cout << "Fail to store Physics Tables for "
474 << particle->GetParticleName()
475 << " and process " << GetProcessName()
476 << " in the directory <" << directory
477 << "> " << G4endl;
478 }
479 }
480 return yes;
481}
482
483//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
484
485G4bool G4VEmProcess::RetrievePhysicsTable(const G4ParticleDefinition* part,
486 const G4String& directory,
487 G4bool ascii)
488{
489 if(1 < verboseLevel) {
490 G4cout << "G4VEmProcess::RetrievePhysicsTable() for "
491 << part->GetParticleName() << " and process "
492 << GetProcessName() << G4endl;
493 }
494 G4bool yes = true;
495
496 if(!buildLambdaTable || particle != part) return yes;
497
498 const G4String particleName = part->GetParticleName();
499 G4String filename;
500
501 filename = GetPhysicsTableFileName(part,directory,"Lambda",ascii);
502 yes = G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTable,
503 filename,ascii);
504 if ( yes ) {
505 if (0 < verboseLevel) {
506 G4cout << "Lambda table for " << particleName << " is Retrieved from <"
507 << filename << ">"
508 << G4endl;
509 }
510 if((G4LossTableManager::Instance())->SplineFlag()) {
511 size_t n = theLambdaTable->length();
512 for(size_t i=0; i<n; i++) {(* theLambdaTable)[i]->SetSpline(true);}
513 }
514 } else {
515 if (1 < verboseLevel) {
516 G4cout << "Lambda table for " << particleName << " in file <"
517 << filename << "> is not exist"
518 << G4endl;
519 }
520 }
521
522 return yes;
523}
524
525//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
526
527void G4VEmProcess::FindLambdaMax()
528{
529 if(1 < verboseLevel) {
530 G4cout << "### G4VEmProcess::FindLambdaMax: "
531 << particle->GetParticleName()
532 << " and process " << GetProcessName() << G4endl;
533 }
534 size_t n = theLambdaTable->length();
535 G4PhysicsVector* pv = (*theLambdaTable)[0];
536 G4double e, s, emax, smax;
537 theEnergyOfCrossSectionMax = new G4double [n];
538 theCrossSectionMax = new G4double [n];
539 G4bool b;
540
541 for (size_t i=0; i<n; i++) {
542 pv = (*theLambdaTable)[i];
543 emax = DBL_MAX;
544 smax = 0.0;
545 if(pv) {
546 size_t nb = pv->GetVectorLength();
547 emax = pv->GetLowEdgeEnergy(nb);
548 smax = 0.0;
549 for (size_t j=0; j<nb; j++) {
550 e = pv->GetLowEdgeEnergy(j);
551 s = pv->GetValue(e,b);
552 if(s > smax) {
553 smax = s;
554 emax = e;
555 }
556 }
557 }
558 theEnergyOfCrossSectionMax[i] = emax;
559 theCrossSectionMax[i] = smax;
560 if(2 < verboseLevel) {
561 G4cout << "For " << particle->GetParticleName()
562 << " Max CS at i= " << i << " emax(MeV)= " << emax/MeV
563 << " lambda= " << smax << G4endl;
564 }
565 }
566}
567
568//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
569
570G4PhysicsVector* G4VEmProcess::LambdaPhysicsVector(const G4MaterialCutsCouple*)
571{
572 G4PhysicsVector* v =
573 new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nLambdaBins);
574 v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
575 return v;
576}
577
578//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
Note: See TracBrowser for help on using the repository browser.