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

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

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