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

Last change on this file since 1228 was 1228, checked in by garnier, 14 years ago

update geant4.9.3 tag
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 $
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.