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

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

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