source: trunk/source/processes/electromagnetic/utils/src/G4VMultipleScattering.cc@ 1315

Last change on this file since 1315 was 1315, checked in by garnier, 15 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 17.3 KB
RevLine 
[819]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//
[1315]26// $Id: G4VMultipleScattering.cc,v 1.82 2010/04/12 11:45:03 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
[819]28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name: G4VMultipleScattering
35//
36// Author: Vladimir Ivanchenko on base of Laszlo Urban code
37//
38// Creation date: 25.03.2003
39//
40// Modifications:
41//
42// 13.04.03 Change printout (V.Ivanchenko)
43// 04-06-03 Fix compilation warnings (V.Ivanchenko)
44// 16-07-03 Use G4VMscModel interface (V.Ivanchenko)
45// 03-11-03 Fix initialisation problem in RetrievePhysicsTable (V.Ivanchenko)
46// 04-11-03 Update PrintInfoDefinition (V.Ivanchenko)
47// 01-03-04 SampleCosineTheta signature changed
48// 22-04-04 SampleCosineTheta signature changed back to original
49// 27-08-04 Add InitialiseForRun method (V.Ivanchneko)
50// 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivantchenko)
51// 11-03-05 Shift verbose level by 1 (V.Ivantchenko)
52// 15-04-05 optimize internal interface (V.Ivanchenko)
53// 15-04-05 remove boundary flag (V.Ivanchenko)
54// 27-10-05 introduce virtual function MscStepLimitation() (V.Ivanchenko)
55// 12-04-07 Add verbosity at destruction (V.Ivanchenko)
56// 27-10-07 Virtual functions moved to source (V.Ivanchenko)
[961]57// 11-03-08 Set skin value does not effect step limit type (V.Ivanchenko)
[1196]58// 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko)
[819]59//
60// Class Description:
61//
62// It is the generic process of multiple scattering it includes common
63// part of calculations for all charged particles
64
65// -------------------------------------------------------------------
66//
67//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
68//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
69
70#include "G4VMultipleScattering.hh"
71#include "G4LossTableManager.hh"
72#include "G4Step.hh"
73#include "G4ParticleDefinition.hh"
74#include "G4VEmFluctuationModel.hh"
75#include "G4DataVector.hh"
76#include "G4PhysicsTable.hh"
77#include "G4PhysicsVector.hh"
78#include "G4PhysicsLogVector.hh"
79#include "G4UnitsTable.hh"
80#include "G4ProductionCutsTable.hh"
81#include "G4Region.hh"
82#include "G4RegionStore.hh"
83#include "G4PhysicsTableHelper.hh"
84#include "G4GenericIon.hh"
85#include "G4Electron.hh"
[1055]86#include "G4EmConfigurator.hh"
[819]87
88//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
89
[961]90G4VMultipleScattering::G4VMultipleScattering(const G4String& name,
91 G4ProcessType type):
92 G4VContinuousDiscreteProcess(name, type),
93 buildLambdaTable(true),
[819]94 theLambdaTable(0),
95 firstParticle(0),
96 stepLimit(fUseSafety),
[1315]97 skin(1.0),
[1055]98 facrange(0.04),
[819]99 facgeom(2.5),
100 latDisplasment(true),
[1196]101 isIon(false),
[961]102 currentParticle(0),
103 currentCouple(0)
[819]104{
[961]105 SetVerboseLevel(1);
106 SetProcessSubType(fMultipleScattering);
107
108 // Size of tables assuming spline
[819]109 minKinEnergy = 0.1*keV;
[1055]110 maxKinEnergy = 10.0*TeV;
111 nBins = 77;
[819]112
[961]113 // default limit on polar angle
114 polarAngleLimit = 0.0;
115
[819]116 pParticleChange = &fParticleChange;
117
118 modelManager = new G4EmModelManager();
119 (G4LossTableManager::Instance())->Register(this);
120
121}
122
123//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
124
125G4VMultipleScattering::~G4VMultipleScattering()
126{
[961]127 if(1 < verboseLevel) {
[819]128 G4cout << "G4VMultipleScattering destruct " << GetProcessName()
129 << G4endl;
[961]130 }
[819]131 delete modelManager;
132 if (theLambdaTable) {
133 theLambdaTable->clearAndDestroy();
134 delete theLambdaTable;
135 }
136 (G4LossTableManager::Instance())->DeRegister(this);
137}
138
139//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
140
[1055]141void G4VMultipleScattering::AddEmModel(G4int order, G4VEmModel* p,
142 const G4Region* region)
143{
144 G4VEmFluctuationModel* fm = 0;
145 modelManager->AddEmModel(order, p, fm, region);
[1315]146 if(p) { p->SetParticleChange(pParticleChange); }
[1055]147}
[1196]148//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
[1055]149
[1196]150void G4VMultipleScattering::SetModel(G4VMscModel* p, G4int index)
151{
152 G4int n = mscModels.size();
153 if(index >= n) { for(G4int i=n; i<=index; ++i) {mscModels.push_back(0);} }
154 mscModels[index] = p;
155}
156
[1055]157//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
158
[1196]159G4VMscModel* G4VMultipleScattering::Model(G4int index)
[1055]160{
[1196]161 G4VMscModel* p = 0;
162 if(index >= 0 && index < G4int(mscModels.size())) { p = mscModels[index]; }
163 return p;
164}
165
166//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
167
168G4VEmModel*
169G4VMultipleScattering::GetModelByIndex(G4int idx, G4bool ver) const
170{
[1055]171 return modelManager->GetModel(idx, ver);
172}
173
174//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
175
[819]176void G4VMultipleScattering::PreparePhysicsTable(const G4ParticleDefinition& part)
177{
178 if (!firstParticle) {
179 currentCouple = 0;
180 if(part.GetParticleType() == "nucleus" &&
[961]181 part.GetParticleSubType() == "generic") {
182 firstParticle = G4GenericIon::GenericIon();
[1196]183 isIon = true;
[961]184 } else {
185 firstParticle = &part;
[1196]186 if(part.GetParticleType() == "nucleus" ||
187 part.GetPDGMass() > GeV) {isIon = true;}
[961]188 }
[1196]189 // limitations for ions
190 if(isIon) {
191 SetStepLimitType(fMinimal);
192 SetLateralDisplasmentFlag(false);
193 SetBuildLambdaTable(false);
194 }
[819]195 currentParticle = &part;
196 }
197
[1315]198 (G4LossTableManager::Instance())->PreparePhysicsTable(&part, this);
199
[819]200 if(1 < verboseLevel) {
201 G4cout << "### G4VMultipleScattering::PrepearPhysicsTable() for "
202 << GetProcessName()
203 << " and particle " << part.GetParticleName()
204 << " local particle " << firstParticle->GetParticleName()
205 << G4endl;
206 }
207
208 if(firstParticle == &part) {
209
210 InitialiseProcess(firstParticle);
211
[1196]212 // initialisation of models
213 G4int nmod = modelManager->NumberOfModels();
214 for(G4int i=0; i<nmod; ++i) {
215 G4VMscModel* msc = static_cast<G4VMscModel*>(modelManager->GetModel(i));
216 if(isIon) {
217 msc->SetStepLimitType(fMinimal);
218 msc->SetLateralDisplasmentFlag(false);
219 msc->SetRangeFactor(0.2);
220 } else {
221 msc->SetStepLimitType(StepLimitType());
222 msc->SetLateralDisplasmentFlag(LateralDisplasmentFlag());
223 msc->SetSkin(Skin());
224 msc->SetRangeFactor(RangeFactor());
225 msc->SetGeomFactor(GeomFactor());
226 }
227 msc->SetPolarAngleLimit(polarAngleLimit);
228 if(msc->HighEnergyLimit() > maxKinEnergy) {
229 msc->SetHighEnergyLimit(maxKinEnergy);
230 }
231 }
[819]232
[1196]233 modelManager->Initialise(firstParticle, G4Electron::Electron(),
234 10.0, verboseLevel);
235
236 // prepare tables
237 if(buildLambdaTable) {
238 theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable);
239 }
[819]240 }
241}
242
243//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
244
[1315]245void G4VMultipleScattering::BuildPhysicsTable(const G4ParticleDefinition& part)
246{
247 G4String num = part.GetParticleName();
248 if(1 < verboseLevel) {
249 G4cout << "### G4VMultipleScattering::BuildPhysicsTable() for "
250 << GetProcessName()
251 << " and particle " << num
252 << G4endl;
253 }
254
255 (G4LossTableManager::Instance())->BuildPhysicsTable(firstParticle);
256
257 if (buildLambdaTable && firstParticle == &part) {
258
259 const G4ProductionCutsTable* theCoupleTable=
260 G4ProductionCutsTable::GetProductionCutsTable();
261 size_t numOfCouples = theCoupleTable->GetTableSize();
262
263 G4bool splineFlag = (G4LossTableManager::Instance())->SplineFlag();
264
265 G4PhysicsLogVector* aVector = 0;
266 G4PhysicsLogVector* bVector = 0;
267
268 for (size_t i=0; i<numOfCouples; ++i) {
269
270 if (theLambdaTable->GetFlag(i)) {
271 // create physics vector and fill it
272 const G4MaterialCutsCouple* couple =
273 theCoupleTable->GetMaterialCutsCouple(i);
274 if(!bVector) {
275 aVector = static_cast<G4PhysicsLogVector*>(PhysicsVector(couple));
276 bVector = aVector;
277 } else {
278 aVector = new G4PhysicsLogVector(*bVector);
279 }
280 //G4PhysicsVector* aVector = PhysicsVector(couple);
281 aVector->SetSpline(splineFlag);
282 modelManager->FillLambdaVector(aVector, couple, false);
283 if(splineFlag) aVector->FillSecondDerivatives();
284 G4PhysicsTableHelper::SetPhysicsVector(theLambdaTable, i, aVector);
285 }
286 }
287
288 if(1 < verboseLevel) {
289 G4cout << "Lambda table is built for "
290 << num
291 << G4endl;
292 }
293 }
294 if(verboseLevel>0 && ( num == "e-" || num == "mu+" ||
295 num == "proton" || num == "pi-" ||
296 num == "GenericIon")) {
297 PrintInfoDefinition();
298 if(2 < verboseLevel && theLambdaTable) G4cout << *theLambdaTable << G4endl;
299 }
300
301 if(1 < verboseLevel) {
302 G4cout << "### G4VMultipleScattering::BuildPhysicsTable() done for "
303 << GetProcessName()
304 << " and particle " << num
305 << G4endl;
306 }
307}
308
309//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
310
[819]311void G4VMultipleScattering::PrintInfoDefinition()
312{
313 if (0 < verboseLevel) {
[961]314 G4cout << G4endl << GetProcessName()
315 << ": for " << firstParticle->GetParticleName()
316 << " SubType= " << GetProcessSubType()
[819]317 << G4endl;
318 if (theLambdaTable) {
319 G4cout << " Lambda tables from "
320 << G4BestUnit(MinKinEnergy(),"Energy")
321 << " to "
322 << G4BestUnit(MaxKinEnergy(),"Energy")
[961]323 << " in " << nBins << " bins, spline: "
324 << (G4LossTableManager::Instance())->SplineFlag()
[819]325 << G4endl;
326 }
327 PrintInfo();
[961]328 modelManager->DumpModelList(verboseLevel);
[819]329 if (2 < verboseLevel) {
330 G4cout << "LambdaTable address= " << theLambdaTable << G4endl;
331 if(theLambdaTable) G4cout << (*theLambdaTable) << G4endl;
332 }
333 }
334}
335
336//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
337
338G4double G4VMultipleScattering::AlongStepGetPhysicalInteractionLength(
339 const G4Track& track,
[961]340 G4double,
[819]341 G4double currentMinimalStep,
[1196]342 G4double&,
[819]343 G4GPILSelection* selection)
344{
345 // get Step limit proposed by the process
[1196]346 *selection = NotCandidateForSelection;
347 G4double x = currentMinimalStep;
348 DefineMaterial(track.GetMaterialCutsCouple());
349 G4double ekin = track.GetKineticEnergy();
350 if(isIon) { ekin *= proton_mass_c2/track.GetDefinition()->GetPDGMass(); }
351 currentModel = static_cast<G4VMscModel*>(SelectModel(ekin));
352 if(x > 0.0 && ekin > 0.0 && currentModel->IsActive(ekin)) {
353 G4double tPathLength =
354 currentModel->ComputeTruePathLengthLimit(track, theLambdaTable, x);
355 if (tPathLength < x) *selection = CandidateForSelection;
356 x = currentModel->ComputeGeomPathLength(tPathLength);
357 // G4cout << "tPathLength= " << tPathLength
358 // << " stepLimit= " << x
359 // << " currentMinimalStep= " << currentMinimalStep<< G4endl;
360 }
361 return x;
[819]362}
363
364//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
365
[1315]366G4double
367G4VMultipleScattering::PostStepGetPhysicalInteractionLength(
368 const G4Track&, G4double, G4ForceCondition* condition)
369{
370 *condition = Forced;
371 return DBL_MAX;
372}
373
374//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
375
376G4VParticleChange*
377G4VMultipleScattering::AlongStepDoIt(const G4Track& track, const G4Step& step)
378{
379 if(currentModel->IsActive(track.GetKineticEnergy())) {
380 fParticleChange.ProposeTrueStepLength(currentModel->ComputeTrueStepLength(step.GetStepLength()));
381 }
382 return &fParticleChange;
383}
384
385//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
386
387G4VParticleChange*
388G4VMultipleScattering::PostStepDoIt(const G4Track& track, const G4Step& step)
389{
390 fParticleChange.Initialize(track);
391 if(currentModel->IsActive(track.GetKineticEnergy())) {
392 currentModel->SampleScattering(track.GetDynamicParticle(),
393 step.GetPostStepPoint()->GetSafety());
394 }
395 return &fParticleChange;
396}
397
398//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
399
[819]400G4double G4VMultipleScattering::GetContinuousStepLimit(
401 const G4Track& track,
402 G4double previousStepSize,
403 G4double currentMinimalStep,
404 G4double& currentSafety)
405{
[1196]406 G4GPILSelection* selection = 0;
407 return AlongStepGetPhysicalInteractionLength(track,previousStepSize,currentMinimalStep,
408 currentSafety, selection);
[819]409}
410
411//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
412
[1315]413G4double G4VMultipleScattering::ContinuousStepLimit(
414 const G4Track& track,
415 G4double previousStepSize,
416 G4double currentMinimalStep,
417 G4double& currentSafety)
418{
419 return GetContinuousStepLimit(track,previousStepSize,currentMinimalStep,
420 currentSafety);
421}
422
423//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
424
[819]425G4double G4VMultipleScattering::GetMeanFreePath(
426 const G4Track&, G4double, G4ForceCondition* condition)
427{
428 *condition = Forced;
429 return DBL_MAX;
430}
431
432//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
433
434G4PhysicsVector* G4VMultipleScattering::PhysicsVector(const G4MaterialCutsCouple* couple)
435{
436 G4int nbins = 3;
437 if( couple->IsUsed() ) nbins = nBins;
438 G4PhysicsVector* v = new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nbins);
[961]439 v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
[819]440 return v;
441}
442
443//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
444
445G4bool G4VMultipleScattering::StorePhysicsTable(const G4ParticleDefinition* part,
446 const G4String& directory,
447 G4bool ascii)
448{
449 G4bool yes = true;
450 if ( theLambdaTable && part == firstParticle) {
451 const G4String name = GetPhysicsTableFileName(part,directory,"Lambda",ascii);
452 G4bool yes = theLambdaTable->StorePhysicsTable(name,ascii);
453
454 if ( yes ) {
455 if ( verboseLevel>0 ) {
456 G4cout << "Physics table are stored for " << part->GetParticleName()
457 << " and process " << GetProcessName()
458 << " in the directory <" << directory
459 << "> " << G4endl;
460 }
461 } else {
462 G4cout << "Fail to store Physics Table for " << part->GetParticleName()
463 << " and process " << GetProcessName()
464 << " in the directory <" << directory
465 << "> " << G4endl;
466 }
467 }
468 return yes;
469}
470
471//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
472
[1196]473G4bool
474G4VMultipleScattering::RetrievePhysicsTable(const G4ParticleDefinition* part,
475 const G4String& directory,
476 G4bool ascii)
[819]477{
478 if(0 < verboseLevel) {
479 G4cout << "G4VMultipleScattering::RetrievePhysicsTable() for "
480 << part->GetParticleName() << " and process "
481 << GetProcessName() << G4endl;
482 }
483 G4bool yes = true;
484
485 if(!buildLambdaTable || firstParticle != part) return yes;
486
487 const G4String particleName = part->GetParticleName();
488
489 G4String filename = GetPhysicsTableFileName(part,directory,"Lambda",ascii);
[1196]490 yes =
491 G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTable,filename,ascii);
[819]492 if ( yes ) {
493 if (0 < verboseLevel) {
[1196]494 G4cout << "Lambda table for " << part->GetParticleName()
495 << " is retrieved from <"
[819]496 << filename << ">"
497 << G4endl;
498 }
[961]499 if((G4LossTableManager::Instance())->SplineFlag()) {
500 size_t n = theLambdaTable->length();
[1196]501 for(size_t i=0; i<n; ++i) {
502 if((* theLambdaTable)[i]) {
503 (* theLambdaTable)[i]->SetSpline(true);
504 }
505 }
[961]506 }
[819]507 } else {
508 if (1 < verboseLevel) {
[1196]509 G4cout << "Lambda table for " << part->GetParticleName()
510 << " in file <"
[819]511 << filename << "> is not exist"
512 << G4endl;
513 }
514 }
515 return yes;
516}
517
518//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
519
Note: See TracBrowser for help on using the repository browser.