source: trunk/source/processes/electromagnetic/utils/src/G4VEmModel.cc

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

update ti head
File size: 9.8 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: G4VEmModel.cc,v 1.37 2010/10/14 16:27:35 vnivanch Exp $
27// GEANT4 tag $Name: emutils-V09-03-23 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name:     G4VEmModel
35//
36// Author:        Vladimir Ivanchenko
37//
38// Creation date: 25.07.2005
39//
40// Modifications:
41// 25.10.2005 Set default highLimit=100.TeV (V.Ivanchenko)
42// 06.02.2006 add method ComputeMeanFreePath() (mma)
43// 16.02.2009 Move implementations of virtual methods to source (VI)
44//
45//
46// Class Description:
47//
48// Abstract interface to energy loss models
49
50// -------------------------------------------------------------------
51//
52
53#include "G4VEmModel.hh"
54#include "G4LossTableManager.hh"
55#include "G4ProductionCutsTable.hh"
56#include "G4ParticleChangeForLoss.hh"
57#include "G4ParticleChangeForGamma.hh"
58
59//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
60//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61
62G4VEmModel::G4VEmModel(const G4String& nam):
63  flucModel(0),anglModel(0), name(nam), lowLimit(0.1*keV), highLimit(100.0*TeV), 
64  eMinActive(0.0),eMaxActive(DBL_MAX),
65  polarAngleLimit(0.0),secondaryThreshold(DBL_MAX),theLPMflag(false),
66  pParticleChange(0),/*nuclearStopping(false),*/
67  currentCouple(0),currentElement(0),
68  nsec(5),flagDeexcitation(false) 
69{
70  xsec.resize(nsec);
71  nSelectors = 0;
72  G4LossTableManager::Instance()->Register(this);
73}
74
75//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
76
77G4VEmModel::~G4VEmModel()
78{
79  G4LossTableManager::Instance()->DeRegister(this);
80  G4int n = elmSelectors.size();
81  if(n > 0) {
82    for(G4int i=0; i<n; ++i) { 
83      delete elmSelectors[i]; 
84    }
85  }
86  delete anglModel;
87}
88
89//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90
91G4ParticleChangeForLoss* G4VEmModel::GetParticleChangeForLoss()
92{
93  G4ParticleChangeForLoss* p = 0;
94  if (pParticleChange) {
95    p = static_cast<G4ParticleChangeForLoss*>(pParticleChange);
96  } else {
97    p = new G4ParticleChangeForLoss();
98    pParticleChange = p;
99  }
100  return p;
101}
102
103//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
104
105G4ParticleChangeForGamma* G4VEmModel::GetParticleChangeForGamma()
106{
107  G4ParticleChangeForGamma* p = 0;
108  if (pParticleChange) {
109    p = static_cast<G4ParticleChangeForGamma*>(pParticleChange);
110  } else {
111    p = new G4ParticleChangeForGamma();
112    pParticleChange = p;
113  }
114  return p;
115}
116
117//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
118
119void G4VEmModel::InitialiseElementSelectors(const G4ParticleDefinition* p, 
120                                            const G4DataVector& cuts)
121{
122  // initialise before run
123  flagDeexcitation = false;
124  G4LossTableManager* man = G4LossTableManager::Instance();
125  G4bool spline = man->SplineFlag();
126
127  // two times less bins because probability functon is normalized
128  // so correspondingly is more smooth
129  G4int nbins = (man->GetNumberOfBinsPerDecade()/3)*
130    G4int(std::log10(highLimit/lowLimit) + 0.5);
131  if(nbins < 5) { nbins = 5; }
132
133  G4ProductionCutsTable* theCoupleTable=
134    G4ProductionCutsTable::GetProductionCutsTable();
135  G4int numOfCouples = theCoupleTable->GetTableSize();
136
137  // prepare vector
138  if(numOfCouples > nSelectors) { 
139    elmSelectors.reserve(numOfCouples); 
140    for(G4int i=nSelectors; i<numOfCouples; ++i) { elmSelectors.push_back(0); }
141    nSelectors = numOfCouples;
142  }
143
144  // initialise vector
145  for(G4int i=0; i<numOfCouples; ++i) {
146    currentCouple = theCoupleTable->GetMaterialCutsCouple(i);
147    const G4Material* material = currentCouple->GetMaterial();
148    G4int idx = currentCouple->GetIndex();
149
150    // selector already exist check if should be deleted
151    G4bool create = true;
152    if(elmSelectors[i]) {
153      if(material == elmSelectors[i]->GetMaterial()) { create = false; }
154      else { delete elmSelectors[i]; }
155    }
156    if(create) {
157      elmSelectors[i] = new G4EmElementSelector(this,material,nbins,
158                                                lowLimit,highLimit,spline);
159    }
160    elmSelectors[i]->Initialise(p, cuts[idx]);
161    //elmSelectors[i]->Dump(p);
162  } 
163}
164
165//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
166
167G4double G4VEmModel::ComputeDEDXPerVolume(const G4Material*,
168                                          const G4ParticleDefinition*,
169                                          G4double,G4double)
170{
171  return 0.0;
172}
173
174//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
175
176G4double G4VEmModel::CrossSectionPerVolume(const G4Material* material,
177                                           const G4ParticleDefinition* p,
178                                           G4double ekin,
179                                           G4double emin,
180                                           G4double emax)
181{
182  SetupForMaterial(p, material, ekin);
183  G4double cross = 0.0;
184  const G4ElementVector* theElementVector = material->GetElementVector();
185  const G4double* theAtomNumDensityVector = material->GetVecNbOfAtomsPerVolume();
186  G4int nelm = material->GetNumberOfElements(); 
187  if(nelm > nsec) {
188    xsec.resize(nelm);
189    nsec = nelm;
190  }
191  for (G4int i=0; i<nelm; i++) {
192    cross += theAtomNumDensityVector[i]*
193      ComputeCrossSectionPerAtom(p,(*theElementVector)[i],ekin,emin,emax);
194    xsec[i] = cross;
195  }
196  return cross;
197}
198
199//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
200
201const G4Element* G4VEmModel::SelectRandomAtom(const G4Material* material,
202                                              const G4ParticleDefinition* pd,
203                                              G4double kinEnergy,
204                                              G4double tcut,
205                                              G4double tmax)
206{
207  const G4ElementVector* theElementVector = material->GetElementVector();
208  G4int n = material->GetNumberOfElements() - 1;
209  currentElement = (*theElementVector)[n];
210  if (n > 0) {
211    G4double x = G4UniformRand()*
212                 G4VEmModel::CrossSectionPerVolume(material,pd,kinEnergy,tcut,tmax);
213    for(G4int i=0; i<n; ++i) {
214      if (x <= xsec[i]) {
215        currentElement = (*theElementVector)[i];
216        break;
217      }
218    }
219  }
220  return currentElement;
221}
222
223//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
224
225G4double G4VEmModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
226                                                G4double, G4double, G4double,
227                                                G4double, G4double)
228{
229  return 0.0;
230}
231
232//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
233
234void G4VEmModel::DefineForRegion(const G4Region*) 
235{}
236
237//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
238
239G4double G4VEmModel::MinEnergyCut(const G4ParticleDefinition*,
240                                  const G4MaterialCutsCouple*)
241{
242  return 0.0;
243}
244
245//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
246
247G4double G4VEmModel::ChargeSquareRatio(const G4Track& track)
248{
249  return GetChargeSquareRatio(track.GetParticleDefinition(), 
250                              track.GetMaterial(), track.GetKineticEnergy());
251}
252
253//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
254
255G4double G4VEmModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
256                                          const G4Material*, G4double)
257{
258  G4double q = p->GetPDGCharge()/CLHEP::eplus;
259  return q*q;
260}
261
262//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
263
264G4double G4VEmModel::GetParticleCharge(const G4ParticleDefinition* p,
265                                       const G4Material*, G4double)
266{
267  return p->GetPDGCharge();
268}
269
270//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
271
272void G4VEmModel::CorrectionsAlongStep(const G4MaterialCutsCouple*,
273                                      const G4DynamicParticle*,
274                                      G4double&,G4double&,G4double)
275{}
276
277//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
278
279void G4VEmModel::SampleDeexcitationAlongStep(const G4Material*,
280                                             const G4Track&,
281                                             G4double& )
282{}
283
284//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
285
286G4double G4VEmModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
287                                        G4double kineticEnergy)
288{
289  return kineticEnergy;
290}
291
292//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
293
294void G4VEmModel::SetupForMaterial(const G4ParticleDefinition*,
295                                  const G4Material*, G4double)
296{}
297
298//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
299
300void 
301G4VEmModel::SetParticleChange(G4VParticleChange* p, G4VEmFluctuationModel* f)
302{
303  if(p && pParticleChange != p) { pParticleChange = p; }
304  flucModel = f;
305}
306
307//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.