| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // -------------------------------------------
28 //
29 // File name: G4LowEnergyBremsstrahlung
30 //
31 // Author: Alessandra Forti, Vladimir I
32 //
33 // Creation date: March 1999
34 //
35 // Modifications:
36 // 18.04.2000 V.L.
37 // - First implementation of continuous energ
38 // 17.02.2000 Veronique Lefebure
39 // - correct bug : the gamma energy was not d
40 // not produced when its energy was < cutFo
41 //
42 // Added Livermore data table construction met
43 // Modified BuildMeanFreePath to read new data
44 // Modified PostStepDoIt to insert sampling wi
45 // Added SelectRandomAtom A. Forti
46 // Added map of the elements A. Forti
47 // 20.09.00 update printout V.Ivanchenko
48 // 24.04.01 V.Ivanchenko remove RogueWave
49 // 29.09.2001 V.Ivanchenko: major revision bas
50 // 10.10.2001 MGP Revision to improve code qua
51 // 18.10.2001 MGP Revision to improve code qua
52 // 28.10.2001 VI Update printout
53 // 29.11.2001 VI New parametrisation
54 // 30.07.2002 VI Fix in restricted energy los
55 // 21.01.2003 VI Cut per region
56 // 21.02.2003 V.Ivanchenko Energy bins for
57 // 28.02.03 V.Ivanchenko Filename is define
58 // 24.03.2003 P.Rodrigues Changes to accommoda
59 // 20.05.2003 MGP Removed memory leak related
60 // 06.11.2003 MGP Improved user interface to
61 //
62 // -------------------------------------------
63
64 #include "G4LowEnergyBremsstrahlung.hh"
65 #include "G4RDeBremsstrahlungSpectrum.hh"
66 #include "G4RDBremsstrahlungCrossSectionHandle
67 #include "G4RDVBremAngularDistribution.hh"
68 #include "G4RDModifiedTsai.hh"
69 #include "G4RDGenerator2BS.hh"
70 #include "G4RDGenerator2BN.hh"
71 #include "G4RDVDataSetAlgorithm.hh"
72 #include "G4RDLogLogInterpolation.hh"
73 #include "G4RDVEMDataSet.hh"
74 #include "G4EnergyLossTables.hh"
75 #include "G4PhysicalConstants.hh"
76 #include "G4SystemOfUnits.hh"
77 #include "G4UnitsTable.hh"
78 #include "G4Electron.hh"
79 #include "G4Gamma.hh"
80 #include "G4ProductionCutsTable.hh"
81
82
83 G4LowEnergyBremsstrahlung::G4LowEnergyBremsstr
84 : G4eLowEnergyLoss(nam),
85 crossSectionHandler(0),
86 theMeanFreePath(0),
87 energySpectrum(0)
88 {
89 cutForPhotons = 0.;
90 verboseLevel = 0;
91 generatorName = "tsai";
92 angularDistribution = new G4RDModifiedTsai("
93 // angularDistribution->PrintGeneratorInforma
94 TsaiAngularDistribution = new G4RDModifiedTs
95 }
96
97 /*
98 G4LowEnergyBremsstrahlung::G4LowEnergyBremsstr
99 : G4eLowEnergyLoss(nam),
100 crossSectionHandler(0),
101 theMeanFreePath(0),
102 energySpectrum(0),
103 angularDistribution(distribution)
104 {
105 cutForPhotons = 0.;
106 verboseLevel = 0;
107
108 angularDistribution->PrintGeneratorInformati
109
110 TsaiAngularDistribution = new G4RDModifiedTs
111 }
112 */
113
114 G4LowEnergyBremsstrahlung::~G4LowEnergyBremsst
115 {
116 if(crossSectionHandler) delete crossSectionH
117 if(energySpectrum) delete energySpectrum;
118 if(theMeanFreePath) delete theMeanFreePath;
119 delete angularDistribution;
120 delete TsaiAngularDistribution;
121 energyBins.clear();
122 }
123
124
125 void G4LowEnergyBremsstrahlung::BuildPhysicsTa
126 {
127 if(verboseLevel > 0) {
128 G4cout << "G4LowEnergyBremsstrahlung::Buil
129 << G4endl;
130 }
131
132 cutForSecondaryPhotons.clear();
133
134 // Create and fill BremsstrahlungParameters
135 if( energySpectrum != 0 ) delete energySpect
136 energyBins.clear();
137 for(size_t i=0; i<15; i++) {
138 G4double x = 0.1*((G4double)i);
139 if(i == 0) x = 0.01;
140 if(i == 10) x = 0.95;
141 if(i == 11) x = 0.97;
142 if(i == 12) x = 0.99;
143 if(i == 13) x = 0.995;
144 if(i == 14) x = 1.0;
145 energyBins.push_back(x);
146 }
147 const G4String dataName("/brem/br-sp.dat");
148 energySpectrum = new G4RDeBremsstrahlungSpec
149
150 if(verboseLevel > 0) {
151 G4cout << "G4LowEnergyBremsstrahlungSpectr
152 << G4endl;
153 }
154
155 // Create and fill G4RDCrossSectionHandler o
156
157 if( crossSectionHandler != 0 ) delete crossS
158 G4RDVDataSetAlgorithm* interpolation = new G
159 G4double lowKineticEnergy = GetLowerBoundEl
160 G4double highKineticEnergy = GetUpperBoundEl
161 G4int totBin = GetNbinEloss();
162 crossSectionHandler = new G4RDBremsstrahlung
163 crossSectionHandler->Initialise(0,lowKinetic
164 crossSectionHandler->LoadShellData("brem/br-
165
166 if (verboseLevel > 0) {
167 G4cout << GetProcessName()
168 << " is created; Cross section data
169 << G4endl;
170 crossSectionHandler->PrintData();
171 G4cout << "Parameters: "
172 << G4endl;
173 energySpectrum->PrintData();
174 }
175
176 // Build loss table for Bremsstrahlung
177
178 BuildLossTable(aParticleType);
179
180 if(verboseLevel > 0) {
181 G4cout << "The loss table is built"
182 << G4endl;
183 }
184
185 if (&aParticleType==G4Electron::Electron())
186
187 RecorderOfElectronProcess[CounterOfElectro
188 CounterOfElectronProcess++;
189 PrintInfoDefinition();
190
191 } else {
192
193 RecorderOfPositronProcess[CounterOfPositro
194 CounterOfPositronProcess++;
195 }
196
197 // Build mean free path data using cut value
198
199 if( theMeanFreePath != 0 ) delete theMeanFre
200 theMeanFreePath = crossSectionHandler->
201 BuildMeanFreePathForMateri
202
203 if(verboseLevel > 0) {
204 G4cout << "The MeanFreePath table is built
205 << G4endl;
206 }
207
208 // Build common DEDX table for all ionisatio
209
210 BuildDEDXTable(aParticleType);
211
212 if(verboseLevel > 0) {
213 G4cout << "G4LowEnergyBremsstrahlung::Buil
214 << G4endl;
215 }
216
217 }
218
219
220 void G4LowEnergyBremsstrahlung::BuildLossTable
221 {
222 // Build table for energy loss due to soft b
223 // the tables are built for *MATERIALS* binn
224
225 G4double lowKineticEnergy = GetLowerBoundEl
226 G4double highKineticEnergy = GetUpperBoundEl
227 size_t totBin = GetNbinEloss();
228
229 // create table
230
231 if (theLossTable) {
232 theLossTable->clearAndDestroy();
233 delete theLossTable;
234 }
235 const G4ProductionCutsTable* theCoupleTable=
236 G4ProductionCutsTable::GetProductionCu
237 size_t numOfCouples = theCoupleTable->GetTab
238 theLossTable = new G4PhysicsTable(numOfCoupl
239
240 // Clean up the vector of cuts
241 cutForSecondaryPhotons.clear();
242
243 // Loop for materials
244
245 for (size_t j=0; j<numOfCouples; j++) {
246
247 // create physics vector and fill it
248 G4PhysicsLogVector* aVector = new G4Physic
249 highKineticEnergy,
250 totBin);
251
252 // get material parameters needed for the
253 const G4MaterialCutsCouple* couple = theCo
254 const G4Material* material= couple->GetMat
255
256 // the cut cannot be below lowest limit
257 G4double tCut = (*(theCoupleTable->GetEner
258 tCut = std::min(highKineticEnergy, tCut);
259 cutForSecondaryPhotons.push_back(tCut);
260
261 const G4ElementVector* theElementVector =
262 size_t NumberOfElements = material->GetNum
263 const G4double* theAtomicNumDensityVector
264 material->GetAtomicNumDensityVector();
265 if(verboseLevel > 1) {
266 G4cout << "Energy loss for material # "
267 << " tCut(keV)= " << tCut/keV
268 << G4endl;
269 }
270
271 // now comes the loop for the kinetic ener
272 for (size_t i = 0; i<totBin; i++) {
273
274 G4double lowEdgeEnergy = aVector->GetLow
275 G4double ionloss = 0.;
276
277 // loop for elements in the material
278 for (size_t iel=0; iel<NumberOfElements;
279 G4int Z = (G4int)((*theElementVector)[
280 G4double e = energySpectrum->AverageEn
281 G4double cs= crossSectionHandler->Find
282 ionloss += e * cs * theAtomicNumDen
283 if(verboseLevel > 1) {
284 G4cout << "Z= " << Z
285 << "; tCut(keV)= " << tCut/ke
286 << "; E(keV)= " << lowEdgeEne
287 << "; Eav(keV)= " << e/keV
288 << "; cs= " << cs
289 << "; loss= " << ionloss
290 << G4endl;
291 }
292 }
293 aVector->PutValue(i,ionloss);
294 }
295 theLossTable->insert(aVector);
296 }
297 }
298
299
300 G4VParticleChange* G4LowEnergyBremsstrahlung::
301 const G4Step& step)
302 {
303 aParticleChange.Initialize(track);
304
305 const G4MaterialCutsCouple* couple = track.G
306 G4double kineticEnergy = track.GetKineticEne
307 G4int index = couple->GetIndex();
308 G4double tCut = cutForSecondaryPhotons[index
309
310 // Control limits
311 if(tCut >= kineticEnergy)
312 return G4VContinuousDiscreteProcess::Post
313
314 G4int Z = crossSectionHandler->SelectRandomA
315
316 G4double tGamma = energySpectrum->SampleEner
317 G4double totalEnergy = kineticEnergy + elect
318 G4double finalEnergy = kineticEnergy - tGamm
319 G4double theta = 0;
320
321 if((kineticEnergy < 1*MeV && kineticEnergy >
322 theta = angularDistribution->PolarAngle(
323 }else{
324 theta = TsaiAngularDistribution->PolarAn
325 }
326
327 G4double phi = twopi * G4UniformRand();
328 G4double dirZ = std::cos(theta);
329 G4double sinTheta = std::sqrt(1. - dirZ*dir
330 G4double dirX = sinTheta*std::cos(phi);
331 G4double dirY = sinTheta*std::sin(phi);
332
333 G4ThreeVector gammaDirection (dirX, dirY, di
334 G4ThreeVector electronDirection = track.GetM
335
336 //
337 // Update the incident particle
338 //
339 gammaDirection.rotateUz(electronDirection);
340
341 // Kinematic problem
342 if (finalEnergy < 0.) {
343 tGamma += finalEnergy;
344 finalEnergy = 0.0;
345 }
346
347 G4double momentum = std::sqrt((totalEnergy +
348
349 G4double finalX = momentum*electronDirection
350 G4double finalY = momentum*electronDirection
351 G4double finalZ = momentum*electronDirection
352
353 aParticleChange.SetNumberOfSecondaries(1);
354 G4double norm = 1./std::sqrt(finalX*finalX +
355 aParticleChange.ProposeMomentumDirection(fin
356 aParticleChange.ProposeEnergy( finalEnergy )
357
358 // create G4DynamicParticle object for the g
359 G4DynamicParticle* aGamma= new G4DynamicPart
360 gammaDirection, tGamma);
361 aParticleChange.AddSecondary(aGamma);
362
363 return G4VContinuousDiscreteProcess::PostSte
364 }
365
366
367 void G4LowEnergyBremsstrahlung::PrintInfoDefin
368 {
369 G4String comments = "Total cross sections fr
370 comments += "\n Gamma energy sampled fr
371 comments += "\n Implementation of the c
372 comments += "\n At present it can be us
373 comments += "in the energy range [250eV,100G
374 comments += "\n The process must work w
375
376 G4cout << G4endl << GetProcessName() << ":
377 }
378
379 G4bool G4LowEnergyBremsstrahlung::IsApplicable
380 {
381 return ( (&particle == G4Electron::Electron
382 }
383
384
385 G4double G4LowEnergyBremsstrahlung::GetMeanFre
386 G4double,
387 G4ForceCondition* cond)
388 {
389 *cond = NotForced;
390 G4int index = (track.GetMaterialCutsCouple()
391 const G4RDVEMDataSet* data = theMeanFreePath
392 G4double meanFreePath = data->FindValue(trac
393 return meanFreePath;
394 }
395
396 void G4LowEnergyBremsstrahlung::SetCutForLowEn
397 {
398 cutForPhotons = cut;
399 }
400
401 void G4LowEnergyBremsstrahlung::SetAngularGene
402 {
403 angularDistribution = distribution;
404 angularDistribution->PrintGeneratorInformati
405 }
406
407 void G4LowEnergyBremsstrahlung::SetAngularGene
408 {
409 if (name == "tsai")
410 {
411 delete angularDistribution;
412 angularDistribution = new G4RDModifiedTs
413 generatorName = name;
414 }
415 else if (name == "2bn")
416 {
417 delete angularDistribution;
418 angularDistribution = new G4RDGenerator2
419 generatorName = name;
420 }
421 else if (name == "2bs")
422 {
423 delete angularDistribution;
424 angularDistribution = new G4RDGenerator
425 generatorName = name;
426 }
427 else
428 {
429 G4Exception("G4LowEnergyBremsstrahlung::
430 "InvalidSetup", FatalExcepti
431 }
432
433 angularDistribution->PrintGeneratorInformati
434 }
435
436