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

Last change on this file since 1199 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 15.7 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: G4VMultipleScattering.cc,v 1.77 2009/10/29 18:07:08 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-03-cand-01 $
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)
57// 11-03-08 Set skin value does not effect step limit type (V.Ivanchenko)
58// 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko)
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"
86#include "G4EmConfigurator.hh"
87
88//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
89
90G4VMultipleScattering::G4VMultipleScattering(const G4String& name,
91 G4ProcessType type):
92 G4VContinuousDiscreteProcess(name, type),
93 buildLambdaTable(true),
94 theLambdaTable(0),
95 firstParticle(0),
96 stepLimit(fUseSafety),
97 skin(3.0),
98 facrange(0.04),
99 facgeom(2.5),
100 latDisplasment(true),
101 isIon(false),
102 currentParticle(0),
103 currentCouple(0)
104{
105 SetVerboseLevel(1);
106 SetProcessSubType(fMultipleScattering);
107
108 // Size of tables assuming spline
109 minKinEnergy = 0.1*keV;
110 maxKinEnergy = 10.0*TeV;
111 nBins = 77;
112
113 // default limit on polar angle
114 polarAngleLimit = 0.0;
115
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{
127 if(1 < verboseLevel) {
128 G4cout << "G4VMultipleScattering destruct " << GetProcessName()
129 << G4endl;
130 }
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
141void G4VMultipleScattering::AddEmModel(G4int order, G4VEmModel* p,
142 const G4Region* region)
143{
144 G4VEmFluctuationModel* fm = 0;
145 modelManager->AddEmModel(order, p, fm, region);
146 if(p) p->SetParticleChange(pParticleChange);
147}
148//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
149
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
157//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
158
159G4VMscModel* G4VMultipleScattering::Model(G4int index)
160{
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{
171 return modelManager->GetModel(idx, ver);
172}
173
174//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
175
176void G4VMultipleScattering::BuildPhysicsTable(const G4ParticleDefinition& part)
177{
178 G4String num = part.GetParticleName();
179 if(1 < verboseLevel) {
180 G4cout << "### G4VMultipleScattering::BuildPhysicsTable() for "
181 << GetProcessName()
182 << " and particle " << num
183 << G4endl;
184 }
185
186 if (buildLambdaTable && firstParticle == &part) {
187
188 const G4ProductionCutsTable* theCoupleTable=
189 G4ProductionCutsTable::GetProductionCutsTable();
190 size_t numOfCouples = theCoupleTable->GetTableSize();
191
192 G4bool splineFlag = (G4LossTableManager::Instance())->SplineFlag();
193
194 G4PhysicsLogVector* aVector = 0;
195 G4PhysicsLogVector* bVector = 0;
196
197 for (size_t i=0; i<numOfCouples; ++i) {
198
199 if (theLambdaTable->GetFlag(i)) {
200 // create physics vector and fill it
201 const G4MaterialCutsCouple* couple =
202 theCoupleTable->GetMaterialCutsCouple(i);
203 if(!bVector) {
204 aVector = static_cast<G4PhysicsLogVector*>(PhysicsVector(couple));
205 bVector = aVector;
206 } else {
207 aVector = new G4PhysicsLogVector(*bVector);
208 }
209 //G4PhysicsVector* aVector = PhysicsVector(couple);
210 aVector->SetSpline(splineFlag);
211 modelManager->FillLambdaVector(aVector, couple, false);
212 if(splineFlag) aVector->FillSecondDerivatives();
213 G4PhysicsTableHelper::SetPhysicsVector(theLambdaTable, i, aVector);
214 }
215 }
216
217 if(1 < verboseLevel) {
218 G4cout << "Lambda table is built for "
219 << num
220 << G4endl;
221 }
222 }
223 if(verboseLevel>0 && ( num == "e-" || num == "mu+" ||
224 num == "proton" || num == "pi-" ||
225 num == "GenericIon")) {
226 PrintInfoDefinition();
227 if(2 < verboseLevel && theLambdaTable) G4cout << *theLambdaTable << G4endl;
228 }
229
230 if(1 < verboseLevel) {
231 G4cout << "### G4VMultipleScattering::BuildPhysicsTable() done for "
232 << GetProcessName()
233 << " and particle " << num
234 << G4endl;
235 }
236}
237
238//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
239
240void G4VMultipleScattering::PreparePhysicsTable(const G4ParticleDefinition& part)
241{
242 if (!firstParticle) {
243 currentCouple = 0;
244 if(part.GetParticleType() == "nucleus" &&
245 part.GetParticleSubType() == "generic") {
246 firstParticle = G4GenericIon::GenericIon();
247 isIon = true;
248 } else {
249 firstParticle = &part;
250 if(part.GetParticleType() == "nucleus" ||
251 part.GetPDGMass() > GeV) {isIon = true;}
252 }
253 // limitations for ions
254 if(isIon) {
255 SetStepLimitType(fMinimal);
256 SetLateralDisplasmentFlag(false);
257 SetBuildLambdaTable(false);
258 }
259 currentParticle = &part;
260 }
261
262 if(1 < verboseLevel) {
263 G4cout << "### G4VMultipleScattering::PrepearPhysicsTable() for "
264 << GetProcessName()
265 << " and particle " << part.GetParticleName()
266 << " local particle " << firstParticle->GetParticleName()
267 << G4endl;
268 }
269
270 (G4LossTableManager::Instance())->EmConfigurator()->AddModels();
271
272 if(firstParticle == &part) {
273
274 InitialiseProcess(firstParticle);
275
276 // initialisation of models
277 G4int nmod = modelManager->NumberOfModels();
278 for(G4int i=0; i<nmod; ++i) {
279 G4VMscModel* msc = static_cast<G4VMscModel*>(modelManager->GetModel(i));
280 if(isIon) {
281 msc->SetStepLimitType(fMinimal);
282 msc->SetLateralDisplasmentFlag(false);
283 msc->SetRangeFactor(0.2);
284 } else {
285 msc->SetStepLimitType(StepLimitType());
286 msc->SetLateralDisplasmentFlag(LateralDisplasmentFlag());
287 msc->SetSkin(Skin());
288 msc->SetRangeFactor(RangeFactor());
289 msc->SetGeomFactor(GeomFactor());
290 }
291 msc->SetPolarAngleLimit(polarAngleLimit);
292 if(msc->HighEnergyLimit() > maxKinEnergy) {
293 msc->SetHighEnergyLimit(maxKinEnergy);
294 }
295 }
296
297 modelManager->Initialise(firstParticle, G4Electron::Electron(),
298 10.0, verboseLevel);
299
300 // prepare tables
301 if(buildLambdaTable) {
302 theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable);
303 }
304 }
305}
306
307//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
308
309void G4VMultipleScattering::PrintInfoDefinition()
310{
311 if (0 < verboseLevel) {
312 G4cout << G4endl << GetProcessName()
313 << ": for " << firstParticle->GetParticleName()
314 << " SubType= " << GetProcessSubType()
315 << G4endl;
316 if (theLambdaTable) {
317 G4cout << " Lambda tables from "
318 << G4BestUnit(MinKinEnergy(),"Energy")
319 << " to "
320 << G4BestUnit(MaxKinEnergy(),"Energy")
321 << " in " << nBins << " bins, spline: "
322 << (G4LossTableManager::Instance())->SplineFlag()
323 << G4endl;
324 }
325 PrintInfo();
326 modelManager->DumpModelList(verboseLevel);
327 if (2 < verboseLevel) {
328 G4cout << "LambdaTable address= " << theLambdaTable << G4endl;
329 if(theLambdaTable) G4cout << (*theLambdaTable) << G4endl;
330 }
331 }
332}
333
334//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
335
336G4double G4VMultipleScattering::AlongStepGetPhysicalInteractionLength(
337 const G4Track& track,
338 G4double,
339 G4double currentMinimalStep,
340 G4double&,
341 G4GPILSelection* selection)
342{
343 // get Step limit proposed by the process
344 *selection = NotCandidateForSelection;
345 G4double x = currentMinimalStep;
346 DefineMaterial(track.GetMaterialCutsCouple());
347 G4double ekin = track.GetKineticEnergy();
348 if(isIon) { ekin *= proton_mass_c2/track.GetDefinition()->GetPDGMass(); }
349 currentModel = static_cast<G4VMscModel*>(SelectModel(ekin));
350 if(x > 0.0 && ekin > 0.0 && currentModel->IsActive(ekin)) {
351 G4double tPathLength =
352 currentModel->ComputeTruePathLengthLimit(track, theLambdaTable, x);
353 if (tPathLength < x) *selection = CandidateForSelection;
354 x = currentModel->ComputeGeomPathLength(tPathLength);
355 // G4cout << "tPathLength= " << tPathLength
356 // << " stepLimit= " << x
357 // << " currentMinimalStep= " << currentMinimalStep<< G4endl;
358 }
359 return x;
360}
361
362//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
363
364G4double G4VMultipleScattering::GetContinuousStepLimit(
365 const G4Track& track,
366 G4double previousStepSize,
367 G4double currentMinimalStep,
368 G4double& currentSafety)
369{
370 G4GPILSelection* selection = 0;
371 return AlongStepGetPhysicalInteractionLength(track,previousStepSize,currentMinimalStep,
372 currentSafety, selection);
373}
374
375//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
376
377G4double G4VMultipleScattering::GetMeanFreePath(
378 const G4Track&, G4double, G4ForceCondition* condition)
379{
380 *condition = Forced;
381 return DBL_MAX;
382}
383
384//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
385
386G4PhysicsVector* G4VMultipleScattering::PhysicsVector(const G4MaterialCutsCouple* couple)
387{
388 G4int nbins = 3;
389 if( couple->IsUsed() ) nbins = nBins;
390 G4PhysicsVector* v = new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nbins);
391 v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
392 return v;
393}
394
395//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
396
397G4bool G4VMultipleScattering::StorePhysicsTable(const G4ParticleDefinition* part,
398 const G4String& directory,
399 G4bool ascii)
400{
401 G4bool yes = true;
402 if ( theLambdaTable && part == firstParticle) {
403 const G4String name = GetPhysicsTableFileName(part,directory,"Lambda",ascii);
404 G4bool yes = theLambdaTable->StorePhysicsTable(name,ascii);
405
406 if ( yes ) {
407 if ( verboseLevel>0 ) {
408 G4cout << "Physics table are stored for " << part->GetParticleName()
409 << " and process " << GetProcessName()
410 << " in the directory <" << directory
411 << "> " << G4endl;
412 }
413 } else {
414 G4cout << "Fail to store Physics Table for " << part->GetParticleName()
415 << " and process " << GetProcessName()
416 << " in the directory <" << directory
417 << "> " << G4endl;
418 }
419 }
420 return yes;
421}
422
423//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
424
425G4bool
426G4VMultipleScattering::RetrievePhysicsTable(const G4ParticleDefinition* part,
427 const G4String& directory,
428 G4bool ascii)
429{
430 if(0 < verboseLevel) {
431 G4cout << "G4VMultipleScattering::RetrievePhysicsTable() for "
432 << part->GetParticleName() << " and process "
433 << GetProcessName() << G4endl;
434 }
435 G4bool yes = true;
436
437 if(!buildLambdaTable || firstParticle != part) return yes;
438
439 const G4String particleName = part->GetParticleName();
440
441 G4String filename = GetPhysicsTableFileName(part,directory,"Lambda",ascii);
442 yes =
443 G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTable,filename,ascii);
444 if ( yes ) {
445 if (0 < verboseLevel) {
446 G4cout << "Lambda table for " << part->GetParticleName()
447 << " is retrieved from <"
448 << filename << ">"
449 << G4endl;
450 }
451 if((G4LossTableManager::Instance())->SplineFlag()) {
452 size_t n = theLambdaTable->length();
453 for(size_t i=0; i<n; ++i) {
454 if((* theLambdaTable)[i]) {
455 (* theLambdaTable)[i]->SetSpline(true);
456 }
457 }
458 }
459 } else {
460 if (1 < verboseLevel) {
461 G4cout << "Lambda table for " << part->GetParticleName()
462 << " in file <"
463 << filename << "> is not exist"
464 << G4endl;
465 }
466 }
467 return yes;
468}
469
470//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
471
Note: See TracBrowser for help on using the repository browser.