Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4VEmModel.cc
Go to the documentation of this file.
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$
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4VEmModel
34 //
35 // Author: Vladimir Ivanchenko
36 //
37 // Creation date: 25.07.2005
38 //
39 // Modifications:
40 // 25.10.2005 Set default highLimit=100.TeV (V.Ivanchenko)
41 // 06.02.2006 add method ComputeMeanFreePath() (mma)
42 // 16.02.2009 Move implementations of virtual methods to source (VI)
43 //
44 //
45 // Class Description:
46 //
47 // Abstract interface to energy loss models
48 
49 // -------------------------------------------------------------------
50 //
51 
52 #include "G4VEmModel.hh"
53 #include "G4LossTableManager.hh"
54 #include "G4ProductionCutsTable.hh"
55 #include "G4ParticleChangeForLoss.hh"
56 #include "G4ParticleChangeForGamma.hh"
57 
58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
60 
61 G4VEmModel::G4VEmModel(const G4String& nam):
62  flucModel(0),anglModel(0), name(nam), lowLimit(0.1*CLHEP::keV),
63  highLimit(100.0*CLHEP::TeV),eMinActive(0.0),eMaxActive(DBL_MAX),
64  polarAngleLimit(CLHEP::pi),secondaryThreshold(DBL_MAX),
65  theLPMflag(false),flagDeexcitation(false),flagForceBuildTable(false),
66  pParticleChange(0),xSectionTable(0),theDensityFactor(0),theDensityIdx(0),
67  fCurrentCouple(0),fCurrentElement(0),
68  nsec(5)
69 {
70  xsec.resize(nsec);
71  nSelectors = 0;
72  G4LossTableManager::Instance()->Register(this);
73 }
74 
75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
76 
77 G4VEmModel::~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  if(xSectionTable) {
88  xSectionTable->clearAndDestroy();
89  delete xSectionTable;
90  }
91 }
92 
93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
94 
95 G4ParticleChangeForLoss* G4VEmModel::GetParticleChangeForLoss()
96 {
97  G4ParticleChangeForLoss* p = 0;
98  if (pParticleChange) {
99  p = static_cast<G4ParticleChangeForLoss*>(pParticleChange);
100  } else {
101  p = new G4ParticleChangeForLoss();
102  pParticleChange = p;
103  }
104  return p;
105 }
106 
107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
108 
109 G4ParticleChangeForGamma* G4VEmModel::GetParticleChangeForGamma()
110 {
111  G4ParticleChangeForGamma* p = 0;
112  if (pParticleChange) {
113  p = static_cast<G4ParticleChangeForGamma*>(pParticleChange);
114  } else {
115  p = new G4ParticleChangeForGamma();
116  pParticleChange = p;
117  }
118  return p;
119 }
120 
121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
122 
123 void G4VEmModel::InitialiseElementSelectors(const G4ParticleDefinition* p,
124  const G4DataVector& cuts)
125 {
126  // initialise before run
127  G4LossTableManager* man = G4LossTableManager::Instance();
128 
129  // using spline for element selectors should be investigated in details
130  // because small number of points may provide biased results
131  // large number of points requires significant increase of memory
132  //G4bool spline = man->SplineFlag();
133  G4bool spline = false;
134 
135  //G4cout << "IES: for " << GetName() << " Emin(MeV)= " << lowLimit/MeV
136  // << " Emax(MeV)= " << highLimit/MeV << G4endl;
137 
138  // two times less bins because probability functon is normalized
139  // so correspondingly is more smooth
140  G4int nbins = G4int(man->GetNumberOfBinsPerDecade()
141  * std::log10(highLimit/lowLimit) / 6.0);
142  if(nbins < 5) { nbins = 5; }
143 
144  G4ProductionCutsTable* theCoupleTable=
145  G4ProductionCutsTable::GetProductionCutsTable();
146  G4int numOfCouples = theCoupleTable->GetTableSize();
147 
148  // prepare vector
149  if(numOfCouples > nSelectors) {
150  elmSelectors.resize(numOfCouples,0);
151  nSelectors = numOfCouples;
152  }
153 
154  // initialise vector
155  for(G4int i=0; i<numOfCouples; ++i) {
156  fCurrentCouple = theCoupleTable->GetMaterialCutsCouple(i);
157  const G4Material* material = fCurrentCouple->GetMaterial();
158  G4int idx = fCurrentCouple->GetIndex();
159 
160  // selector already exist check if should be deleted
161  G4bool create = true;
162  if(elmSelectors[i]) {
163  if(material == elmSelectors[i]->GetMaterial()) { create = false; }
164  else { delete elmSelectors[i]; }
165  }
166  if(create) {
167  elmSelectors[i] = new G4EmElementSelector(this,material,nbins,
168  lowLimit,highLimit,spline);
169  }
170  elmSelectors[i]->Initialise(p, cuts[idx]);
171  //elmSelectors[i]->Dump(p);
172  }
173 }
174 
175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
176 
177 G4double G4VEmModel::ComputeDEDXPerVolume(const G4Material*,
178  const G4ParticleDefinition*,
179  G4double,G4double)
180 {
181  return 0.0;
182 }
183 
184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
185 
186 G4double G4VEmModel::CrossSectionPerVolume(const G4Material* material,
187  const G4ParticleDefinition* p,
188  G4double ekin,
189  G4double emin,
190  G4double emax)
191 {
192  SetupForMaterial(p, material, ekin);
193  G4double cross = 0.0;
194  const G4ElementVector* theElementVector = material->GetElementVector();
195  const G4double* theAtomNumDensityVector = material->GetVecNbOfAtomsPerVolume();
196  G4int nelm = material->GetNumberOfElements();
197  if(nelm > nsec) {
198  xsec.resize(nelm);
199  nsec = nelm;
200  }
201  for (G4int i=0; i<nelm; i++) {
202  cross += theAtomNumDensityVector[i]*
203  ComputeCrossSectionPerAtom(p,(*theElementVector)[i],ekin,emin,emax);
204  xsec[i] = cross;
205  }
206  return cross;
207 }
208 
209 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
210 
211 void G4VEmModel::StartTracking(G4Track*)
212 {}
213 
214 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
215 
216 const G4Element* G4VEmModel::SelectRandomAtom(const G4Material* material,
217  const G4ParticleDefinition* pd,
218  G4double kinEnergy,
219  G4double tcut,
220  G4double tmax)
221 {
222  const G4ElementVector* theElementVector = material->GetElementVector();
223  G4int n = material->GetNumberOfElements() - 1;
224  fCurrentElement = (*theElementVector)[n];
225  if (n > 0) {
226  G4double x = G4UniformRand()*
227  G4VEmModel::CrossSectionPerVolume(material,pd,kinEnergy,tcut,tmax);
228  for(G4int i=0; i<n; ++i) {
229  if (x <= xsec[i]) {
230  fCurrentElement = (*theElementVector)[i];
231  break;
232  }
233  }
234  }
235  return fCurrentElement;
236 }
237 
238 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
239 
240 G4double G4VEmModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
241  G4double, G4double, G4double,
242  G4double, G4double)
243 {
244  return 0.0;
245 }
246 
247 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
248 
249 void G4VEmModel::DefineForRegion(const G4Region*)
250 {}
251 
252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
253 
254 G4double G4VEmModel::ChargeSquareRatio(const G4Track& track)
255 {
256  return GetChargeSquareRatio(track.GetParticleDefinition(),
257  track.GetMaterial(), track.GetKineticEnergy());
258 }
259 
260 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
261 
262 G4double G4VEmModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
263  const G4Material*, G4double)
264 {
265  G4double q = p->GetPDGCharge()/CLHEP::eplus;
266  return q*q;
267 }
268 
269 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
270 
271 G4double G4VEmModel::GetParticleCharge(const G4ParticleDefinition* p,
272  const G4Material*, G4double)
273 {
274  return p->GetPDGCharge();
275 }
276 
277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
278 
279 void G4VEmModel::CorrectionsAlongStep(const G4MaterialCutsCouple*,
280  const G4DynamicParticle*,
281  G4double&,G4double&,G4double)
282 {}
283 
284 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
285 
286 G4double G4VEmModel::Value(const G4MaterialCutsCouple* couple,
287  const G4ParticleDefinition* p, G4double e)
288 {
289  fCurrentCouple = couple;
290  return e*e*CrossSectionPerVolume(couple->GetMaterial(),p,e,0.0,DBL_MAX);
291 }
292 
293 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
294 
295 G4double G4VEmModel::MinPrimaryEnergy(const G4Material*,
296  const G4ParticleDefinition*)
297 {
298  return 0.0;
299 }
300 
301 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
302 
303 G4double G4VEmModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
304  G4double kineticEnergy)
305 {
306  return kineticEnergy;
307 }
308 
309 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
310 
311 void G4VEmModel::SetupForMaterial(const G4ParticleDefinition*,
312  const G4Material*, G4double)
313 {}
314 
315 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
316 
317 void
318 G4VEmModel::SetParticleChange(G4VParticleChange* p, G4VEmFluctuationModel* f)
319 {
320  if(p && pParticleChange != p) { pParticleChange = p; }
321  flucModel = f;
322 }
323 
324 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
325 
326 void G4VEmModel::SetCrossSectionTable(G4PhysicsTable* p)
327 {
328  if(p != xSectionTable) {
329  if(xSectionTable) {
330  xSectionTable->clearAndDestroy();
331  delete xSectionTable;
332  }
333  xSectionTable = p;
334  }
335 }
336 
337 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Generated on Sat May 25 2013 14:33:39 for Geant4 by   doxygen 1.8.4