| 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$
26 // 27 //
27 // ------------------------------------------- << 28 // Author: Luciano Pandola
>> 29 // on base of G4LowEnergyBremsstrahlung developed by A.Forti and V.Ivanchenko
28 // 30 //
29 // GEANT4 Class file << 31 // History:
>> 32 // --------
>> 33 // 03 Mar 2009 L Pandola Migration from process to model
>> 34 // 12 Apr 2009 V Ivanchenko Cleanup initialisation and generation of secondaries:
>> 35 // - apply internal high-energy limit only in constructor
>> 36 // - do not apply low-energy limit (default is 0)
>> 37 // - added MinEnergyCut method
>> 38 // - do not change track status
>> 39 // - do not initialize element selectors
>> 40 // - use cut value from the interface
>> 41 // - fixed bug in sampling of angles between keV and MeV
>> 42 // 19 May 2009 L Pandola Explicitely set to zero pointers deleted in
>> 43 // Initialise(), since they might be checked later on
30 // 44 //
31 // <<
32 // File name: G4LivermoreBremsstrahlungMod <<
33 // <<
34 // Author: Vladimir Ivanchenko use inhe <<
35 // base class implementing ultr <<
36 // model <<
37 // <<
38 // Creation date: 04.10.2011 <<
39 // <<
40 // Modifications: <<
41 // <<
42 // ------------------------------------------- <<
43 // <<
44 //....oooOO0OOooo........oooOO0OOooo........oo <<
45 //....oooOO0OOooo........oooOO0OOooo........oo <<
46 45
47 #include "G4LivermoreBremsstrahlungModel.hh" 46 #include "G4LivermoreBremsstrahlungModel.hh"
48 #include "G4PhysicalConstants.hh" 47 #include "G4PhysicalConstants.hh"
49 #include "G4SystemOfUnits.hh" 48 #include "G4SystemOfUnits.hh"
50 #include "G4Electron.hh" << 49 #include "G4ParticleDefinition.hh"
51 #include "G4Positron.hh" << 50 #include "G4MaterialCutsCouple.hh"
52 #include "G4Gamma.hh" << 51
53 #include "Randomize.hh" << 52 #include "G4DynamicParticle.hh"
54 #include "G4AutoLock.hh" <<
55 #include "G4Material.hh" <<
56 #include "G4Element.hh" 53 #include "G4Element.hh"
57 #include "G4ElementVector.hh" << 54 #include "G4Gamma.hh"
58 #include "G4ProductionCutsTable.hh" << 55 #include "G4Electron.hh"
59 #include "G4ParticleChangeForLoss.hh" << 56 #include "G4SemiLogInterpolation.hh"
>> 57 //
>> 58 #include "G4VEmAngularDistribution.hh"
>> 59 #include "G4ModifiedTsai.hh"
60 #include "G4Generator2BS.hh" 60 #include "G4Generator2BS.hh"
61 << 61 //#include "G4Generator2BN.hh"
62 #include "G4Physics2DVector.hh" << 62 //
63 #include "G4Exp.hh" << 63 #include "G4BremsstrahlungCrossSectionHandler.hh"
64 #include "G4Log.hh" << 64 //
65 << 65 #include "G4VEnergySpectrum.hh"
66 #include "G4ios.hh" << 66 #include "G4eBremsstrahlungSpectrum.hh"
67 #include <fstream> << 67 #include "G4VEMDataSet.hh"
68 #include <iomanip> <<
69 68
70 //....oooOO0OOooo........oooOO0OOooo........oo 69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71 70
72 namespace { G4Mutex LivermoreBremsstrahlungMod <<
73 using namespace std; <<
74 71
75 << 72 G4LivermoreBremsstrahlungModel::G4LivermoreBremsstrahlungModel(const G4ParticleDefinition*,
76 G4Physics2DVector* G4LivermoreBremsstrahlungMo << 73 const G4String& nam)
77 G4double G4LivermoreBremsstrahlungModel::ylimi << 74 :G4VEmModel(nam),fParticleChange(0),isInitialised(false),
78 G4double G4LivermoreBremsstrahlungModel::expnu << 75 crossSectionHandler(0),energySpectrum(0)
79 <<
80 static const G4double emaxlog = 4*G4Log(10.); <<
81 static const G4double alpha = CLHEP::twopi*CLH <<
82 static const G4double epeaklimit= 300*CLHEP::M <<
83 static const G4double elowlimit = 20*CLHEP::ke <<
84 <<
85 G4LivermoreBremsstrahlungModel::G4LivermoreBre <<
86 const G4ParticleDefinition* p, const G4Strin <<
87 : G4eBremsstrahlungRelModel(p,nam),useBicubi <<
88 { 76 {
89 SetLowEnergyLimit(10.0*eV); << 77 fIntrinsicLowEnergyLimit = 10.0*eV;
>> 78 fIntrinsicHighEnergyLimit = 100.0*GeV;
>> 79 fNBinEnergyLoss = 360;
>> 80 // SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
>> 81 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
>> 82 //
>> 83 verboseLevel = 0;
90 SetAngularDistribution(new G4Generator2BS()) 84 SetAngularDistribution(new G4Generator2BS());
>> 85 //
>> 86 //generatorName = "tsai";
>> 87 //angularDistribution = new G4ModifiedTsai("TsaiGenerator"); //default generator
>> 88 //
>> 89 //TsaiAngularDistribution = new G4ModifiedTsai("TsaiGenerator");
>> 90 //
91 } 91 }
92 92
93 //....oooOO0OOooo........oooOO0OOooo........oo 93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
94 94
95 G4LivermoreBremsstrahlungModel::~G4LivermoreBr 95 G4LivermoreBremsstrahlungModel::~G4LivermoreBremsstrahlungModel()
96 { 96 {
97 if(IsMaster()) { << 97 if (crossSectionHandler) delete crossSectionHandler;
98 for(size_t i=0; i<101; ++i) { << 98 if (energySpectrum) delete energySpectrum;
99 if(dataSB[i]) { << 99 energyBins.clear();
100 delete dataSB[i]; << 100 //delete angularDistribution;
101 dataSB[i] = nullptr; << 101 //delete TsaiAngularDistribution;
102 } <<
103 } <<
104 } <<
105 } 102 }
106 103
107 //....oooOO0OOooo........oooOO0OOooo........oo 104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
108 105
109 void G4LivermoreBremsstrahlungModel::Initialis << 106 void G4LivermoreBremsstrahlungModel::Initialise(const G4ParticleDefinition* particle,
110 const G4DataVector& cuts) 107 const G4DataVector& cuts)
111 { 108 {
112 // Access to elements << 109 //Check that the Livermore Bremsstrahlung is NOT attached to e+
113 if(IsMaster()) { << 110 if (particle != G4Electron::Electron())
114 // check environment variable << 111 {
115 // Build the complete string identifying t << 112 G4Exception("G4LivermoreBremsstrahlungModel::Initialise",
116 const char* path = G4FindDataDir("G4LEDATA << 113 "em0002",FatalException,"Livermore Bremsstrahlung Model is applicable only to electrons");
117 << 114 }
118 const G4ElementTable* theElmTable = G4Elem << 115 //Prepare energy spectrum
119 size_t numOfElm = G4Element::GetNumberOfEl << 116 if (energySpectrum)
120 if(numOfElm > 0) { << 117 {
121 for(size_t i=0; i<numOfElm; ++i) { << 118 delete energySpectrum;
122 G4int Z = (*theElmTable)[i]->GetZasInt(); << 119 energySpectrum = 0;
123 if(Z < 1) { Z = 1; } << 120 }
124 else if(Z > 100) { Z = 100; } << 121
125 //G4cout << "Z= " << Z << G4endl; << 122 energyBins.clear();
126 // Initialisation << 123 for(size_t i=0; i<15; i++)
127 if(!dataSB[Z]) { ReadData(Z, path); } << 124 {
128 } << 125 G4double x = 0.1*((G4double)i);
>> 126 if(i == 0) x = 0.01;
>> 127 if(i == 10) x = 0.95;
>> 128 if(i == 11) x = 0.97;
>> 129 if(i == 12) x = 0.99;
>> 130 if(i == 13) x = 0.995;
>> 131 if(i == 14) x = 1.0;
>> 132 energyBins.push_back(x);
>> 133 }
>> 134 const G4String dataName("/brem/br-sp.dat");
>> 135 energySpectrum = new G4eBremsstrahlungSpectrum(energyBins,dataName);
>> 136
>> 137 if (verboseLevel > 0)
>> 138 G4cout << "G4eBremsstrahlungSpectrum is initialized" << G4endl;
>> 139
>> 140 //Initialize cross section handler
>> 141 if (crossSectionHandler)
>> 142 {
>> 143 delete crossSectionHandler;
>> 144 crossSectionHandler = 0;
>> 145 }
>> 146 G4VDataSetAlgorithm* interpolation = 0;//new G4SemiLogInterpolation();
>> 147 crossSectionHandler = new G4BremsstrahlungCrossSectionHandler(energySpectrum,interpolation);
>> 148 crossSectionHandler->Initialise(0,LowEnergyLimit(),HighEnergyLimit(),
>> 149 fNBinEnergyLoss);
>> 150 crossSectionHandler->Clear();
>> 151 crossSectionHandler->LoadShellData("brem/br-cs-");
>> 152 //This is used to retrieve cross section values later on
>> 153 G4VEMDataSet* p = crossSectionHandler->BuildMeanFreePathForMaterials(&cuts);
>> 154 delete p;
>> 155
>> 156 if (verboseLevel > 0)
>> 157 {
>> 158 G4cout << "Livermore Bremsstrahlung model is initialized " << G4endl
>> 159 << "Energy range: "
>> 160 << LowEnergyLimit() / keV << " keV - "
>> 161 << HighEnergyLimit() / GeV << " GeV"
>> 162 << G4endl;
129 } 163 }
130 } << 164
131 G4eBremsstrahlungRelModel::Initialise(p, cut << 165 if (verboseLevel > 1)
>> 166 {
>> 167 G4cout << "Cross section data: " << G4endl;
>> 168 crossSectionHandler->PrintData();
>> 169 G4cout << "Parameters: " << G4endl;
>> 170 energySpectrum->PrintData();
>> 171 }
>> 172
>> 173 if(isInitialised) return;
>> 174 fParticleChange = GetParticleChangeForLoss();
>> 175 isInitialised = true;
132 } 176 }
133 177
134 //....oooOO0OOooo........oooOO0OOooo........oo 178 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
135 179
136 G4String G4LivermoreBremsstrahlungModel::Direc << 180 G4double G4LivermoreBremsstrahlungModel::MinEnergyCut(const G4ParticleDefinition*,
>> 181 const G4MaterialCutsCouple*)
137 { 182 {
138 return "/livermore/brem/br"; << 183 return 250.*eV;
139 } 184 }
140 185
141 //....oooOO0OOooo........oooOO0OOooo........oo 186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
142 187
143 void G4LivermoreBremsstrahlungModel::ReadData( << 188 G4double
>> 189 G4LivermoreBremsstrahlungModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
>> 190 G4double energy,
>> 191 G4double Z, G4double,
>> 192 G4double cutEnergy,
>> 193 G4double)
144 { 194 {
145 if(dataSB[Z]) { return; } << 195 G4int iZ = (G4int) Z;
146 const char* datadir = path; << 196 if (!crossSectionHandler)
147 << 197 {
148 if(nullptr == datadir) { << 198 G4Exception("G4LivermoreBremsstrahlungModel::ComputeCrossSectionPerAtom",
149 datadir = G4FindDataDir("G4LEDATA"); << 199 "em1007",FatalException,"The cross section handler is not correctly initialized");
150 if(!datadir) { << 200 return 0;
151 G4Exception("G4LivermoreBremsstrahlungMo << 201 }
152 FatalException,"Environment variable G4L << 202
153 return; << 203 //The cut is already included in the crossSectionHandler
>> 204 G4double cs =
>> 205 crossSectionHandler->GetCrossSectionAboveThresholdForElement(energy,cutEnergy,iZ);
>> 206
>> 207 if (verboseLevel > 1)
>> 208 {
>> 209 G4cout << "G4LivermoreBremsstrahlungModel " << G4endl;
>> 210 G4cout << "Cross section for gamma emission > " << cutEnergy/keV << " keV at " <<
>> 211 energy/keV << " keV and Z = " << iZ << " --> " << cs/barn << " barn" << G4endl;
154 } 212 }
155 } << 213 return cs;
156 std::ostringstream ost; <<
157 ost << datadir << DirectoryPath() << Z; <<
158 std::ifstream fin(ost.str().c_str()); <<
159 if( !fin.is_open()) { <<
160 G4ExceptionDescription ed; <<
161 ed << "Bremsstrahlung data file <" << ost. <<
162 << "> is not opened!"; <<
163 G4Exception("G4LivermoreBremsstrahlungMode <<
164 FatalException,ed, <<
165 "G4LEDATA version should be G4EMLOW8.0 or <<
166 return; <<
167 } <<
168 //G4cout << "G4LivermoreBremsstrahlungModel <<
169 // << ">" << G4endl; <<
170 G4Physics2DVector* v = new G4Physics2DVector <<
171 if(v->Retrieve(fin)) { <<
172 if(useBicubicInterpolation) { v->SetBicubi <<
173 dataSB[Z] = v; <<
174 ylimit[Z] = v->Value(0.97, emaxlog, idx, i <<
175 } else { <<
176 G4ExceptionDescription ed; <<
177 ed << "Bremsstrahlung data file <" << ost. <<
178 << "> is not retrieved!"; <<
179 G4Exception("G4LivermoreBremsstrahlungMode <<
180 FatalException,ed, <<
181 "G4LEDATA version should be G4EMLOW8.0 or <<
182 delete v; <<
183 } <<
184 } 214 }
185 215
>> 216
186 //....oooOO0OOooo........oooOO0OOooo........oo 217 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
187 218
188 G4double << 219 G4double G4LivermoreBremsstrahlungModel::ComputeDEDXPerVolume(const G4Material* material,
189 G4LivermoreBremsstrahlungModel::ComputeDXSecti << 220 const G4ParticleDefinition* ,
>> 221 G4double kineticEnergy,
>> 222 G4double cutEnergy)
190 { 223 {
191 if(gammaEnergy < 0.0 || fPrimaryKinEnergy <= << 224 G4double sPower = 0.0;
192 G4double x = gammaEnergy/fPrimaryKinEnergy; << 225
193 G4double y = G4Log(fPrimaryKinEnergy/MeV); << 226 const G4ElementVector* theElementVector = material->GetElementVector();
194 G4int Z = fCurrentIZ; << 227 size_t NumberOfElements = material->GetNumberOfElements() ;
195 << 228 const G4double* theAtomicNumDensityVector =
196 //G4cout << "G4LivermoreBremsstrahlungModel: << 229 material->GetAtomicNumDensityVector();
197 // << " x= " << x << " y= " << y << " " << << 230
198 if(!dataSB[Z]) { InitialiseForElement(0, Z); << 231 // loop for elements in the material
199 << 232 for (size_t iel=0; iel<NumberOfElements; iel++ )
200 G4double invb2 = fPrimaryTotalEnergy*fPrimar << 233 {
201 *(fPrimaryKinEnergy + 2.*fP << 234 G4int iZ = (G4int)((*theElementVector)[iel]->GetZ());
202 G4double cross = dataSB[Z]->Value(x,y,idx,id << 235 G4double e = energySpectrum->AverageEnergy(iZ, 0.0,cutEnergy,
203 << 236 kineticEnergy);
204 if(!fIsElectron) { << 237 G4double cs= crossSectionHandler->FindValue(iZ,kineticEnergy);
205 G4double invbeta1 = sqrt(invb2); << 238 sPower += e * cs * theAtomicNumDensityVector[iel];
206 G4double e2 = fPrimaryKinEnergy - gammaEne <<
207 if(e2 > 0.0) { <<
208 G4double invbeta2 = (e2 + fPrimaryPartic <<
209 /sqrt(e2*(e2 + 2.*fP <<
210 G4double xxx = alpha*fCurrentIZ*(invbeta <<
211 if(xxx < expnumlim) { cross = 0.0; } <<
212 else { cross *= G4Exp(xxx); } <<
213 } else { <<
214 cross = 0.0; <<
215 } 239 }
216 } <<
217 240
218 return cross; << 241 if (verboseLevel > 2)
>> 242 {
>> 243 G4cout << "G4LivermoreBremsstrahlungModel " << G4endl;
>> 244 G4cout << "Stopping power < " << cutEnergy/keV << " keV at " <<
>> 245 kineticEnergy/keV << " keV = " << sPower/(keV/mm) << " keV/mm" << G4endl;
>> 246 }
>> 247
>> 248 return sPower;
219 } 249 }
220 250
221 //....oooOO0OOooo........oooOO0OOooo........oo 251 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
222 252
223 void << 253 void G4LivermoreBremsstrahlungModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
224 G4LivermoreBremsstrahlungModel::SampleSecondar << 254 const G4MaterialCutsCouple* couple,
225 std::v << 255 const G4DynamicParticle* aDynamicParticle,
226 const G4MaterialCutsCouple* couple, << 256 G4double energyCut,
227 const G4DynamicParticle* dp, << 257 G4double)
228 G4double cutEnergy, <<
229 G4double maxEnergy) <<
230 { 258 {
231 G4double kineticEnergy = dp->GetKineticEnerg << 259
232 G4double cut = std::min(cutEnergy, kineticE << 260 G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
233 G4double emax = std::min(maxEnergy, kineticE <<
234 if(cut >= emax) { return; } <<
235 // sets total energy, kinetic energy and den <<
236 SetupForMaterial(fPrimaryParticle, couple->G <<
237 <<
238 const G4Element* elm = <<
239 SelectRandomAtom(couple,fPrimaryParticle,k <<
240 fCurrentIZ = elm->GetZasInt(); <<
241 G4int Z = fCurrentIZ; <<
242 <<
243 G4double totMomentum = sqrt(kineticEnergy*(f <<
244 /* <<
245 G4cout << "G4LivermoreBremsstrahlungModel::S <<
246 << kineticEnergy/MeV <<
247 << " Z= " << Z << " cut(MeV)= " << cut/MeV <<
248 << " emax(MeV)= " << emax/MeV << " corr= " <<
249 */ <<
250 G4double xmin = G4Log(cut*cut + fDensityCorr <<
251 G4double xmax = G4Log(emax*emax + fDensityC <<
252 G4double y = G4Log(kineticEnergy/MeV); <<
253 <<
254 G4double gammaEnergy, v; <<
255 <<
256 // majoranta <<
257 G4double x0 = cut/kineticEnergy; <<
258 G4double vmax = dataSB[Z]->Value(x0, y, idx, <<
259 <<
260 // majoranta corrected for e- <<
261 if(fIsElectron && x0 < 0.97 && <<
262 ((kineticEnergy > epeaklimit) || (kinetic <<
263 G4double ylim = std::min(ylimit[Z],1.1*dat <<
264 if(ylim > vmax) { vmax = ylim; } <<
265 } <<
266 if(x0 < 0.05) { vmax *= 1.2; } <<
267 <<
268 do { <<
269 //++ncount; <<
270 G4double x = G4Exp(xmin + G4UniformRand()* <<
271 if(x < 0.0) { x = 0.0; } <<
272 gammaEnergy = sqrt(x); <<
273 G4double x1 = gammaEnergy/kineticEnergy; <<
274 v = dataSB[Z]->Value(x1, y, idx, idy); <<
275 <<
276 // correction for positrons <<
277 if(!fIsElectron) { <<
278 G4double e1 = kineticEnergy - cut; <<
279 G4double invbeta1 = (e1 + fPrimaryPartic <<
280 /sqrt(e1*(e1 + 2*fPr <<
281 G4double e2 = kineticEnergy - gammaEnerg <<
282 G4double invbeta2 = (e2 + fPrimaryPartic <<
283 /sqrt(e2*(e2 + 2*fPr <<
284 G4double xxx = twopi*fine_structure_cons <<
285 <<
286 if(xxx < expnumlim) { v = 0.0; } <<
287 else { v *= G4Exp(xxx); } <<
288 } <<
289 <<
290 if (v > 1.05*vmax && nwarn < 5) { <<
291 ++nwarn; <<
292 G4ExceptionDescription ed; <<
293 ed << "### G4LivermoreBremsstrahlungMode <<
294 << v << " > " << vmax << " by " << v/vmax <<
295 << " Egamma(MeV)= " << gammaEnergy <<
296 << " Ee(MeV)= " << kineticEnergy <<
297 << " Z= " << Z << " " << fPrimaryParticle- <<
298 <<
299 if ( 20 == nwarn ) { <<
300 ed << "\n ### G4LivermoreBremsstrahlungModel <<
301 } <<
302 G4Exception("G4LivermoreBremsstrahlungMo <<
303 JustWarning, ed,""); <<
304 261
>> 262 // this is neede for pathalogical cases of no ionisation
>> 263 if (kineticEnergy <= fIntrinsicLowEnergyLimit)
>> 264 {
>> 265 fParticleChange->SetProposedKineticEnergy(0.);
>> 266 fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy);
>> 267 return;
305 } 268 }
306 } while (v < vmax*G4UniformRand()); <<
307 269
308 // << 270 //Sample material
309 // angles of the emitted gamma. ( Z - axis a << 271 G4int Z = crossSectionHandler->SelectRandomAtom(couple, kineticEnergy);
310 // use general interface <<
311 // <<
312 272
313 G4ThreeVector gammaDirection = << 273 //Sample gamma energy
314 GetAngularDistribution()->SampleDirection( << 274 G4double tGamma = energySpectrum->SampleEnergy(Z, energyCut, kineticEnergy, kineticEnergy);
315 Z, couple->GetMaterial()); << 275 //nothing happens
316 << 276 if (tGamma == 0.) { return; }
317 // create G4DynamicParticle object for the G << 277
318 G4DynamicParticle* gamma = << 278 G4double totalEnergy = kineticEnergy + electron_mass_c2;
319 new G4DynamicParticle(fGammaParticle,gamma << 279 G4double finalEnergy = kineticEnergy - tGamma; // electron final energy
320 vdp->push_back(gamma); << 280
321 << 281 //Sample gamma direction
322 G4ThreeVector direction = (totMomentum*dp->G << 282 G4ThreeVector gammaDirection =
323 - gammaEnergy*gammaDirection).unit( << 283 GetAngularDistribution()->SampleDirection(aDynamicParticle,
324 << 284 totalEnergy-tGamma,
325 /* << 285 Z,
326 G4cout << "### G4SBModel: v= " << 286 couple->GetMaterial());
327 << " Eg(MeV)= " << gammaEnergy << 287
328 << " Ee(MeV)= " << kineticEnergy << 288 G4ThreeVector electronDirection = aDynamicParticle->GetMomentumDirection();
329 << " DirE " << direction << " DirG " << gam << 289
330 << G4endl; << 290 //Update the incident particle
331 */ << 291 if (finalEnergy < 0.)
332 // energy of primary << 292 {
333 G4double finalE = kineticEnergy - gammaEnerg << 293 // Kinematic problem
334 << 294 tGamma = kineticEnergy;
335 // stop tracking and create new secondary in << 295 fParticleChange->SetProposedKineticEnergy(0.);
336 if(gammaEnergy > SecondaryThreshold()) { << 296 }
337 fParticleChange->ProposeTrackStatus(fStopA << 297 else
338 fParticleChange->SetProposedKineticEnergy( << 298 {
339 G4DynamicParticle* el = << 299 G4double momentum = std::sqrt((totalEnergy + electron_mass_c2)*kineticEnergy);
340 new G4DynamicParticle(const_cast<G4Parti << 300 G4double finalX = momentum*electronDirection.x() - tGamma*gammaDirection.x();
341 direction, finalE); << 301 G4double finalY = momentum*electronDirection.y() - tGamma*gammaDirection.y();
342 vdp->push_back(el); << 302 G4double finalZ = momentum*electronDirection.z() - tGamma*gammaDirection.z();
343 << 303 G4double norm = 1./std::sqrt(finalX*finalX + finalY*finalY + finalZ*finalZ);
344 // continue tracking << 304
345 } else { << 305 fParticleChange->ProposeMomentumDirection(finalX*norm, finalY*norm, finalZ*norm);
346 fParticleChange->SetProposedMomentumDirect << 306 fParticleChange->SetProposedKineticEnergy(finalEnergy);
347 fParticleChange->SetProposedKineticEnergy( << 307 }
348 } <<
349 } <<
350 308
351 //....oooOO0OOooo........oooOO0OOooo........oo << 309 //Generate the bremsstrahlung gamma
>> 310 G4DynamicParticle* aGamma= new G4DynamicParticle (G4Gamma::Gamma(),
>> 311 gammaDirection, tGamma);
>> 312 fvect->push_back(aGamma);
>> 313
>> 314 if (verboseLevel > 1)
>> 315 {
>> 316 G4cout << "-----------------------------------------------------------" << G4endl;
>> 317 G4cout << "Energy balance from G4LivermoreBremsstrahlung" << G4endl;
>> 318 G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl;
>> 319 G4cout << "-----------------------------------------------------------" << G4endl;
>> 320 G4cout << "Outgoing primary energy: " << finalEnergy/keV << " keV" << G4endl;
>> 321 G4cout << "Gamma ray " << tGamma/keV << " keV" << G4endl;
>> 322 G4cout << "Total final state: " << (finalEnergy+tGamma)/keV << " keV" << G4endl;
>> 323 G4cout << "-----------------------------------------------------------" << G4endl;
>> 324 }
>> 325 if (verboseLevel > 0)
>> 326 {
>> 327 G4double energyDiff = std::fabs(finalEnergy+tGamma-kineticEnergy);
>> 328 if (energyDiff > 0.05*keV)
>> 329 G4cout << "G4LivermoreBremsstrahlung WARNING: problem with energy conservation: "
>> 330 << (finalEnergy+tGamma)/keV << " keV (final) vs. "
>> 331 << kineticEnergy/keV << " keV (initial)" << G4endl;
>> 332 }
>> 333 return;
>> 334 }
352 335
353 void G4LivermoreBremsstrahlungModel::Initialis << 336 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
354 const G4P << 337 /*
355 G4int Z) << 338 void
>> 339 G4LivermoreBremsstrahlungModel::SetAngularGenerator(G4VBremAngularDistribution* distribution)
356 { 340 {
357 G4AutoLock l(&LivermoreBremsstrahlungModelMu << 341 if(angularDistribution == distribution) return;
358 if(!dataSB[Z]) { ReadData(Z); } << 342 if(angularDistribution) delete angularDistribution;
359 l.unlock(); << 343 angularDistribution = distribution;
>> 344 angularDistribution->PrintGeneratorInformation();
360 } 345 }
>> 346 */
>> 347 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 348 /*
>> 349 void G4LivermoreBremsstrahlungModel::SetAngularGenerator(const G4String& theGenName)
>> 350 {
>> 351 if(theGenName == generatorName) return;
>> 352 if (theGenName == "tsai")
>> 353 {
>> 354 delete angularDistribution;
>> 355 angularDistribution = new G4ModifiedTsai("TsaiGenerator");
>> 356 generatorName = theGenName;
>> 357 }
>> 358 else if (theGenName == "2bn")
>> 359 {
>> 360 delete angularDistribution;
>> 361 angularDistribution = new G4Generator2BN("2BNGenerator");
>> 362 generatorName = theGenName;
>> 363 }
>> 364 else if (theGenName == "2bs")
>> 365 {
>> 366 delete angularDistribution;
>> 367 angularDistribution = new G4Generator2BS("2BSGenerator");
>> 368 generatorName = theGenName;
>> 369 }
>> 370 else
>> 371 {
>> 372 G4cout << "### G4LivermoreBremsstrahlungModel::SetAngularGenerator WARNING:"
>> 373 << " generator <" << theGenName << "> is not known" << G4endl;
>> 374 return;
361 375
362 //....oooOO0OOooo........oooOO0OOooo........oo << 376 }
>> 377
>> 378 angularDistribution->PrintGeneratorInformation();
>> 379 }
>> 380 */
363 381