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

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

maj sur la beta de geant 4.9.3
File size: 14.4 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.66 2009/05/27 11:55:02 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-03-beta-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//
59// Class Description:
60//
61// It is the generic process of multiple scattering it includes common
62// part of calculations for all charged particles
63
64// -------------------------------------------------------------------
65//
66//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
67//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
68
69#include "G4VMultipleScattering.hh"
70#include "G4LossTableManager.hh"
71#include "G4Step.hh"
72#include "G4ParticleDefinition.hh"
73#include "G4VEmFluctuationModel.hh"
74#include "G4DataVector.hh"
75#include "G4PhysicsTable.hh"
76#include "G4PhysicsVector.hh"
77#include "G4PhysicsLogVector.hh"
78#include "G4UnitsTable.hh"
79#include "G4ProductionCutsTable.hh"
80#include "G4Region.hh"
81#include "G4RegionStore.hh"
82#include "G4PhysicsTableHelper.hh"
83#include "G4GenericIon.hh"
84#include "G4Electron.hh"
85#include "G4EmConfigurator.hh"
86
87//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
88
89G4VMultipleScattering::G4VMultipleScattering(const G4String& name,
90 G4ProcessType type):
91 G4VContinuousDiscreteProcess(name, type),
92 buildLambdaTable(true),
93 theLambdaTable(0),
94 firstParticle(0),
95 stepLimit(fUseSafety),
96 skin(3.0),
97 facrange(0.04),
98 facgeom(2.5),
99 latDisplasment(true),
100 currentParticle(0),
101 currentCouple(0)
102{
103 SetVerboseLevel(1);
104 SetProcessSubType(fMultipleScattering);
105
106 // Size of tables assuming spline
107 minKinEnergy = 0.1*keV;
108 maxKinEnergy = 10.0*TeV;
109 nBins = 77;
110
111 // default limit on polar angle
112 polarAngleLimit = 0.0;
113
114 pParticleChange = &fParticleChange;
115
116 modelManager = new G4EmModelManager();
117 (G4LossTableManager::Instance())->Register(this);
118
119}
120
121//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
122
123G4VMultipleScattering::~G4VMultipleScattering()
124{
125 if(1 < verboseLevel) {
126 G4cout << "G4VMultipleScattering destruct " << GetProcessName()
127 << G4endl;
128 }
129 delete modelManager;
130 if (theLambdaTable) {
131 theLambdaTable->clearAndDestroy();
132 delete theLambdaTable;
133 }
134 (G4LossTableManager::Instance())->DeRegister(this);
135}
136
137//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
138
139void G4VMultipleScattering::AddEmModel(G4int order, G4VEmModel* p,
140 const G4Region* region)
141{
142 G4VEmFluctuationModel* fm = 0;
143 modelManager->AddEmModel(order, p, fm, region);
144 if(p) p->SetParticleChange(pParticleChange);
145}
146
147//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
148
149G4VEmModel* G4VMultipleScattering::GetModelByIndex(G4int idx, G4bool ver)
150{
151 return modelManager->GetModel(idx, ver);
152}
153
154//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
155
156void G4VMultipleScattering::BuildPhysicsTable(const G4ParticleDefinition& part)
157{
158 G4String num = part.GetParticleName();
159 if(1 < verboseLevel) {
160 G4cout << "### G4VMultipleScattering::BuildPhysicsTable() for "
161 << GetProcessName()
162 << " and particle " << num
163 << G4endl;
164 }
165
166 if (buildLambdaTable && firstParticle == &part) {
167
168 const G4ProductionCutsTable* theCoupleTable=
169 G4ProductionCutsTable::GetProductionCutsTable();
170 size_t numOfCouples = theCoupleTable->GetTableSize();
171
172 for (size_t i=0; i<numOfCouples; i++) {
173
174 if (theLambdaTable->GetFlag(i)) {
175 // create physics vector and fill it
176 const G4MaterialCutsCouple* couple =
177 theCoupleTable->GetMaterialCutsCouple(i);
178 G4PhysicsVector* aVector = PhysicsVector(couple);
179 modelManager->FillLambdaVector(aVector, couple, false);
180 G4PhysicsTableHelper::SetPhysicsVector(theLambdaTable, i, aVector);
181 }
182 }
183
184 if(1 < verboseLevel) {
185 G4cout << "Lambda table is built for "
186 << num
187 << G4endl;
188 }
189 }
190 if(verboseLevel>0 && ( num == "e-" || num == "mu+" ||
191 num == "proton" || num == "pi-" ||
192 num == "GenericIon")) {
193 PrintInfoDefinition();
194 if(2 < verboseLevel && theLambdaTable) G4cout << *theLambdaTable << G4endl;
195 }
196
197 if(1 < verboseLevel) {
198 G4cout << "### G4VMultipleScattering::BuildPhysicsTable() done for "
199 << GetProcessName()
200 << " and particle " << num
201 << G4endl;
202 }
203}
204
205//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
206
207void G4VMultipleScattering::PreparePhysicsTable(const G4ParticleDefinition& part)
208{
209 if (!firstParticle) {
210 currentCouple = 0;
211 if(part.GetParticleType() == "nucleus" &&
212 part.GetParticleSubType() == "generic") {
213 firstParticle = G4GenericIon::GenericIon();
214 } else {
215 firstParticle = &part;
216 }
217
218 currentParticle = &part;
219 }
220
221 if(1 < verboseLevel) {
222 G4cout << "### G4VMultipleScattering::PrepearPhysicsTable() for "
223 << GetProcessName()
224 << " and particle " << part.GetParticleName()
225 << " local particle " << firstParticle->GetParticleName()
226 << G4endl;
227 }
228
229 (G4LossTableManager::Instance())->EmConfigurator()->AddModels();
230
231 if(firstParticle == &part) {
232
233 InitialiseProcess(firstParticle);
234 if(buildLambdaTable)
235 theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable);
236 const G4DataVector* theCuts =
237 modelManager->Initialise(firstParticle,
238 G4Electron::Electron(),
239 10.0, verboseLevel);
240
241 if(2 < verboseLevel) G4cout << theCuts << G4endl;
242
243 }
244}
245
246//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
247
248void G4VMultipleScattering::PrintInfoDefinition()
249{
250 if (0 < verboseLevel) {
251 G4cout << G4endl << GetProcessName()
252 << ": for " << firstParticle->GetParticleName()
253 << " SubType= " << GetProcessSubType()
254 << G4endl;
255 if (theLambdaTable) {
256 G4cout << " Lambda tables from "
257 << G4BestUnit(MinKinEnergy(),"Energy")
258 << " to "
259 << G4BestUnit(MaxKinEnergy(),"Energy")
260 << " in " << nBins << " bins, spline: "
261 << (G4LossTableManager::Instance())->SplineFlag()
262 << G4endl;
263 }
264 PrintInfo();
265 modelManager->DumpModelList(verboseLevel);
266 if (2 < verboseLevel) {
267 G4cout << "LambdaTable address= " << theLambdaTable << G4endl;
268 if(theLambdaTable) G4cout << (*theLambdaTable) << G4endl;
269 }
270 }
271}
272
273//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
274
275G4double G4VMultipleScattering::AlongStepGetPhysicalInteractionLength(
276 const G4Track& track,
277 G4double,
278 G4double currentMinimalStep,
279 G4double& currentSafety,
280 G4GPILSelection* selection)
281{
282 // get Step limit proposed by the process
283 valueGPILSelectionMSC = NotCandidateForSelection;
284 G4double steplength = GetMscContinuousStepLimit(track,
285 track.GetKineticEnergy(),
286 currentMinimalStep,
287 currentSafety);
288 // G4cout << "StepLimit= " << steplength << G4endl;
289 // set return value for G4GPILSelection
290 *selection = valueGPILSelectionMSC;
291 return steplength;
292}
293
294//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
295
296G4double G4VMultipleScattering::PostStepGetPhysicalInteractionLength(
297 const G4Track&, G4double, G4ForceCondition* condition)
298{
299 *condition = Forced;
300 return DBL_MAX;
301}
302
303//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
304
305G4double G4VMultipleScattering::GetContinuousStepLimit(
306 const G4Track& track,
307 G4double previousStepSize,
308 G4double currentMinimalStep,
309 G4double& currentSafety)
310{
311 return GetMscContinuousStepLimit(track,previousStepSize,currentMinimalStep,
312 currentSafety);
313}
314
315//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
316
317G4double G4VMultipleScattering::GetMeanFreePath(
318 const G4Track&, G4double, G4ForceCondition* condition)
319{
320 *condition = Forced;
321 return DBL_MAX;
322}
323
324//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
325
326G4VParticleChange* G4VMultipleScattering::AlongStepDoIt(const G4Track&,
327 const G4Step& step)
328{
329 fParticleChange.ProposeTrueStepLength(
330 currentModel->ComputeTrueStepLength(step.GetStepLength()));
331 return &fParticleChange;
332}
333
334//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
335
336G4VParticleChange* G4VMultipleScattering::PostStepDoIt(const G4Track& track,
337 const G4Step& step)
338{
339 fParticleChange.Initialize(track);
340 currentModel->SampleScattering(track.GetDynamicParticle(),
341 step.GetPostStepPoint()->GetSafety());
342 return &fParticleChange;
343}
344
345//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
346
347G4PhysicsVector* G4VMultipleScattering::PhysicsVector(const G4MaterialCutsCouple* couple)
348{
349 G4int nbins = 3;
350 if( couple->IsUsed() ) nbins = nBins;
351 G4PhysicsVector* v = new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nbins);
352 v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
353 return v;
354}
355
356//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
357
358G4bool G4VMultipleScattering::StorePhysicsTable(const G4ParticleDefinition* part,
359 const G4String& directory,
360 G4bool ascii)
361{
362 G4bool yes = true;
363 if ( theLambdaTable && part == firstParticle) {
364 const G4String name = GetPhysicsTableFileName(part,directory,"Lambda",ascii);
365 G4bool yes = theLambdaTable->StorePhysicsTable(name,ascii);
366
367 if ( yes ) {
368 if ( verboseLevel>0 ) {
369 G4cout << "Physics table are stored for " << part->GetParticleName()
370 << " and process " << GetProcessName()
371 << " in the directory <" << directory
372 << "> " << G4endl;
373 }
374 } else {
375 G4cout << "Fail to store Physics Table for " << part->GetParticleName()
376 << " and process " << GetProcessName()
377 << " in the directory <" << directory
378 << "> " << G4endl;
379 }
380 }
381 return yes;
382}
383
384//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
385
386G4bool G4VMultipleScattering::RetrievePhysicsTable(const G4ParticleDefinition* part,
387 const G4String& directory,
388 G4bool ascii)
389{
390 if(0 < verboseLevel) {
391// G4cout << "========================================================" << G4endl;
392 G4cout << "G4VMultipleScattering::RetrievePhysicsTable() for "
393 << part->GetParticleName() << " and process "
394 << GetProcessName() << G4endl;
395 }
396 G4bool yes = true;
397
398 if(!buildLambdaTable || firstParticle != part) return yes;
399
400 const G4String particleName = part->GetParticleName();
401
402 G4String filename = GetPhysicsTableFileName(part,directory,"Lambda",ascii);
403 yes = G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTable,filename,ascii);
404 if ( yes ) {
405 if (0 < verboseLevel) {
406 G4cout << "Lambda table for " << part->GetParticleName() << " is retrieved from <"
407 << filename << ">"
408 << G4endl;
409 }
410 if((G4LossTableManager::Instance())->SplineFlag()) {
411 size_t n = theLambdaTable->length();
412 for(size_t i=0; i<n; i++) {(* theLambdaTable)[i]->SetSpline(true);}
413 }
414 } else {
415 if (1 < verboseLevel) {
416 G4cout << "Lambda table for " << part->GetParticleName() << " in file <"
417 << filename << "> is not exist"
418 << G4endl;
419 }
420 }
421 return yes;
422}
423
424//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
425
Note: See TracBrowser for help on using the repository browser.