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

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

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 17.3 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.82 2010/04/12 11:45:03 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-04-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// 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(1.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::PreparePhysicsTable(const G4ParticleDefinition& part)
177{
178  if (!firstParticle) {
179    currentCouple = 0;
180    if(part.GetParticleType() == "nucleus" && 
181       part.GetParticleSubType() == "generic") {
182      firstParticle = G4GenericIon::GenericIon();
183      isIon = true;
184    } else {
185      firstParticle = &part;
186      if(part.GetParticleType() == "nucleus" || 
187         part.GetPDGMass() > GeV) {isIon = true;} 
188    } 
189    // limitations for ions
190    if(isIon) {
191      SetStepLimitType(fMinimal);
192      SetLateralDisplasmentFlag(false);
193      SetBuildLambdaTable(false);
194    }
195    currentParticle = &part;
196  }
197
198  (G4LossTableManager::Instance())->PreparePhysicsTable(&part, this);
199
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
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    }
232
233    modelManager->Initialise(firstParticle, G4Electron::Electron(), 
234                             10.0, verboseLevel);
235
236    // prepare tables
237    if(buildLambdaTable) {
238      theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable);
239    }
240  }
241}
242
243//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
244
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
311void G4VMultipleScattering::PrintInfoDefinition()
312{
313  if (0 < verboseLevel) {
314    G4cout << G4endl << GetProcessName() 
315           << ":   for " << firstParticle->GetParticleName()
316           << "    SubType= " << GetProcessSubType() 
317           << G4endl;
318    if (theLambdaTable) {
319      G4cout << "      Lambda tables from "
320             << G4BestUnit(MinKinEnergy(),"Energy")
321             << " to "
322             << G4BestUnit(MaxKinEnergy(),"Energy")
323             << " in " << nBins << " bins, spline: " 
324             << (G4LossTableManager::Instance())->SplineFlag()
325             << G4endl;
326    }
327    PrintInfo();
328    modelManager->DumpModelList(verboseLevel);
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,
340                             G4double,
341                             G4double currentMinimalStep,
342                             G4double&,
343                             G4GPILSelection* selection)
344{
345  // get Step limit proposed by the process
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;
362}
363
364//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
365
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
400G4double G4VMultipleScattering::GetContinuousStepLimit(
401                                       const G4Track& track,
402                                       G4double previousStepSize,
403                                       G4double currentMinimalStep,
404                                       G4double& currentSafety)
405{
406  G4GPILSelection* selection = 0;
407  return AlongStepGetPhysicalInteractionLength(track,previousStepSize,currentMinimalStep,
408                                               currentSafety, selection);
409}
410
411//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
412
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
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);
439  v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
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
473G4bool
474G4VMultipleScattering::RetrievePhysicsTable(const G4ParticleDefinition* part,
475                                            const G4String& directory,
476                                            G4bool ascii)
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);
490  yes = 
491    G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTable,filename,ascii);
492  if ( yes ) {
493    if (0 < verboseLevel) {
494        G4cout << "Lambda table for " << part->GetParticleName() 
495               << " is retrieved from <"
496               << filename << ">"
497               << G4endl;
498    }
499    if((G4LossTableManager::Instance())->SplineFlag()) {
500      size_t n = theLambdaTable->length();
501      for(size_t i=0; i<n; ++i) {
502        if((* theLambdaTable)[i]) {
503          (* theLambdaTable)[i]->SetSpline(true);
504        }
505      }
506    }
507  } else {
508    if (1 < verboseLevel) {
509        G4cout << "Lambda table for " << part->GetParticleName() 
510               << " in file <"
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.