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