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