| 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: G4TablesForExtrapolator.cc 75145 2013-10-28 18:34:48Z vnivanch $
>> 27 //
26 //-------------------------------------------- 28 //---------------------------------------------------------------------------
27 // 29 //
28 // ClassName: G4TablesForExtrapolator 30 // ClassName: G4TablesForExtrapolator
29 // 31 //
30 // Description: This class provide calculatio 32 // Description: This class provide calculation of energy loss, fluctuation,
31 // and msc angle 33 // and msc angle
32 // 34 //
33 // Author: 09.12.04 V.Ivanchenko 35 // Author: 09.12.04 V.Ivanchenko
34 // 36 //
35 // Modification: 37 // Modification:
36 // 08-04-05 Rename Propogator -> Extrapolator 38 // 08-04-05 Rename Propogator -> Extrapolator (V.Ivanchenko)
37 // 16-03-06 Add muon tables and fix bug in uni 39 // 16-03-06 Add muon tables and fix bug in units (V.Ivanchenko)
38 // 21-03-06 Add verbosity defined in the const 40 // 21-03-06 Add verbosity defined in the constructor and Initialisation
39 // start only when first public metho 41 // start only when first public method is called (V.Ivanchenko)
40 // 03-05-06 Remove unused pointer G4Material* 42 // 03-05-06 Remove unused pointer G4Material* from number of methods (VI)
41 // 12-05-06 SEt linLossLimit=0.001 (VI) 43 // 12-05-06 SEt linLossLimit=0.001 (VI)
42 // 44 //
43 //-------------------------------------------- 45 //----------------------------------------------------------------------------
44 // 46 //
45 47
46 //....oooOO0OOooo........oooOO0OOooo........oo 48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47 49
48 #include "G4TablesForExtrapolator.hh" 50 #include "G4TablesForExtrapolator.hh"
49 #include "G4PhysicalConstants.hh" 51 #include "G4PhysicalConstants.hh"
50 #include "G4SystemOfUnits.hh" 52 #include "G4SystemOfUnits.hh"
51 #include "G4LossTableManager.hh" <<
52 #include "G4PhysicsLogVector.hh" 53 #include "G4PhysicsLogVector.hh"
53 #include "G4ParticleDefinition.hh" 54 #include "G4ParticleDefinition.hh"
54 #include "G4Material.hh" 55 #include "G4Material.hh"
55 #include "G4MaterialCutsCouple.hh" 56 #include "G4MaterialCutsCouple.hh"
56 #include "G4Electron.hh" 57 #include "G4Electron.hh"
57 #include "G4Positron.hh" 58 #include "G4Positron.hh"
58 #include "G4Proton.hh" 59 #include "G4Proton.hh"
59 #include "G4MuonPlus.hh" 60 #include "G4MuonPlus.hh"
60 #include "G4MuonMinus.hh" 61 #include "G4MuonMinus.hh"
61 #include "G4EmParameters.hh" 62 #include "G4EmParameters.hh"
62 #include "G4MollerBhabhaModel.hh" 63 #include "G4MollerBhabhaModel.hh"
63 #include "G4BetheBlochModel.hh" 64 #include "G4BetheBlochModel.hh"
64 #include "G4eBremsstrahlungRelModel.hh" 65 #include "G4eBremsstrahlungRelModel.hh"
65 #include "G4MuPairProductionModel.hh" 66 #include "G4MuPairProductionModel.hh"
66 #include "G4MuBremsstrahlungModel.hh" 67 #include "G4MuBremsstrahlungModel.hh"
67 #include "G4ProductionCuts.hh" 68 #include "G4ProductionCuts.hh"
68 #include "G4LossTableBuilder.hh" 69 #include "G4LossTableBuilder.hh"
69 #include "G4WentzelVIModel.hh" 70 #include "G4WentzelVIModel.hh"
70 #include "G4ios.hh" <<
71 71
72 //....oooOO0OOooo........oooOO0OOooo........oo 72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73 73
74 G4TablesForExtrapolator::G4TablesForExtrapolat 74 G4TablesForExtrapolator::G4TablesForExtrapolator(G4int verb, G4int bins,
75 75 G4double e1, G4double e2)
76 : emin(e1), emax(e2), verbose(verb), nbins(b << 76 : verbose(verb), nbins(bins), emin(e1), emax(e2)
77 { 77 {
78 electron = G4Electron::Electron(); << 78 Initialisation();
79 positron = G4Positron::Positron(); <<
80 proton = G4Proton::Proton(); <<
81 muonPlus = G4MuonPlus::MuonPlus(); <<
82 muonMinus= G4MuonMinus::MuonMinus(); <<
83 } 79 }
84 80
85 //....oooOO0OOooo........oooOO0OOooo........oo 81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
86 82
87 G4TablesForExtrapolator::~G4TablesForExtrapola << 83 G4TablesForExtrapolator:: ~G4TablesForExtrapolator()
88 { 84 {
89 if(nullptr != dedxElectron) { << 85 for(G4int i=0; i<nmat; i++) {delete couples[i];}
90 dedxElectron->clearAndDestroy(); << 86
91 delete dedxElectron; << 87 dedxElectron->clearAndDestroy();
92 } << 88 dedxPositron->clearAndDestroy();
93 if(nullptr != dedxPositron) { << 89 dedxProton->clearAndDestroy();
94 dedxPositron->clearAndDestroy(); << 90 dedxMuon->clearAndDestroy();
95 delete dedxPositron; << 91 rangeElectron->clearAndDestroy();
96 } << 92 rangePositron->clearAndDestroy();
97 if(nullptr != dedxProton) { << 93 rangeProton->clearAndDestroy();
98 dedxProton->clearAndDestroy(); << 94 rangeMuon->clearAndDestroy();
99 delete dedxProton; << 95 invRangeElectron->clearAndDestroy();
100 } << 96 invRangePositron->clearAndDestroy();
101 if(nullptr != dedxMuon) { << 97 invRangeProton->clearAndDestroy();
102 dedxMuon->clearAndDestroy(); << 98 invRangeMuon->clearAndDestroy();
103 delete dedxMuon; << 99 mscElectron->clearAndDestroy();
104 } << 100
105 if(nullptr != rangeElectron) { << 101 delete dedxElectron;
106 rangeElectron->clearAndDestroy(); << 102 delete dedxPositron;
107 delete rangeElectron; << 103 delete dedxProton;
108 } << 104 delete dedxMuon;
109 if(nullptr != rangePositron) { << 105 delete rangeElectron;
110 rangePositron->clearAndDestroy(); << 106 delete rangePositron;
111 delete rangePositron; << 107 delete rangeProton;
112 } << 108 delete rangeMuon;
113 if(nullptr != rangeProton) { << 109 delete invRangeElectron;
114 rangeProton->clearAndDestroy(); << 110 delete invRangePositron;
115 delete rangeProton; << 111 delete invRangeProton;
116 } << 112 delete invRangeMuon;
117 if(nullptr != rangeMuon) { << 113 delete mscElectron;
118 rangeMuon->clearAndDestroy(); <<
119 delete rangeMuon; <<
120 } <<
121 if(nullptr != invRangeElectron) { <<
122 invRangeElectron->clearAndDestroy(); <<
123 delete invRangeElectron; <<
124 } <<
125 if(nullptr != invRangePositron) { <<
126 invRangePositron->clearAndDestroy(); <<
127 delete invRangePositron; <<
128 } <<
129 if(nullptr != invRangeProton) { <<
130 invRangeProton->clearAndDestroy(); <<
131 delete invRangeProton; <<
132 } <<
133 if(nullptr != invRangeMuon) { <<
134 invRangeMuon->clearAndDestroy(); <<
135 delete invRangeMuon; <<
136 } <<
137 if(nullptr != mscElectron) { <<
138 mscElectron->clearAndDestroy(); <<
139 delete mscElectron; <<
140 } <<
141 delete pcuts; 114 delete pcuts;
142 delete builder; <<
143 } 115 }
144 116
145 //....oooOO0OOooo........oooOO0OOooo........oo 117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
146 118
147 const G4PhysicsTable* 119 const G4PhysicsTable*
148 G4TablesForExtrapolator::GetPhysicsTable(ExtTa 120 G4TablesForExtrapolator::GetPhysicsTable(ExtTableType type) const
149 { 121 {
150 const G4PhysicsTable* table = nullptr; 122 const G4PhysicsTable* table = nullptr;
151 switch (type) 123 switch (type)
152 { 124 {
153 case fDedxElectron: 125 case fDedxElectron:
154 table = dedxElectron; 126 table = dedxElectron;
155 break; 127 break;
156 case fDedxPositron: 128 case fDedxPositron:
157 table = dedxPositron; 129 table = dedxPositron;
158 break; 130 break;
159 case fDedxProton: 131 case fDedxProton:
160 table = dedxProton; 132 table = dedxProton;
161 break; 133 break;
162 case fDedxMuon: 134 case fDedxMuon:
163 table = dedxMuon; 135 table = dedxMuon;
164 break; 136 break;
165 case fRangeElectron: 137 case fRangeElectron:
166 table = rangeElectron; 138 table = rangeElectron;
167 break; 139 break;
168 case fRangePositron: 140 case fRangePositron:
169 table = rangePositron; 141 table = rangePositron;
170 break; 142 break;
171 case fRangeProton: 143 case fRangeProton:
172 table = rangeProton; 144 table = rangeProton;
173 break; 145 break;
174 case fRangeMuon: 146 case fRangeMuon:
175 table = rangeMuon; 147 table = rangeMuon;
176 break; 148 break;
177 case fInvRangeElectron: 149 case fInvRangeElectron:
178 table = invRangeElectron; 150 table = invRangeElectron;
179 break; 151 break;
180 case fInvRangePositron: 152 case fInvRangePositron:
181 table = invRangePositron; 153 table = invRangePositron;
182 break; 154 break;
183 case fInvRangeProton: 155 case fInvRangeProton:
184 table = invRangeProton; 156 table = invRangeProton;
185 break; 157 break;
186 case fInvRangeMuon: 158 case fInvRangeMuon:
187 table = invRangeMuon; 159 table = invRangeMuon;
188 break; 160 break;
189 case fMscElectron: 161 case fMscElectron:
190 table = mscElectron; 162 table = mscElectron;
191 } 163 }
192 return table; 164 return table;
193 } 165 }
194 166
195 //....oooOO0OOooo........oooOO0OOooo........oo 167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
196 168
197 void G4TablesForExtrapolator::Initialisation() 169 void G4TablesForExtrapolator::Initialisation()
198 { 170 {
199 if(verbose>1) { 171 if(verbose>1) {
200 G4cout << "### G4TablesForExtrapolator::In 172 G4cout << "### G4TablesForExtrapolator::Initialisation" << G4endl;
201 } 173 }
202 G4int num = (G4int)G4Material::GetNumberOfMa << 174 currentParticle = nullptr;
203 if(nmat == num) { return; } <<
204 nmat = num; <<
205 cuts.resize(nmat, DBL_MAX); <<
206 couples.resize(nmat, nullptr); <<
207 175
>> 176 electron = G4Electron::Electron();
>> 177 positron = G4Positron::Positron();
>> 178 proton = G4Proton::Proton();
>> 179 muonPlus = G4MuonPlus::MuonPlus();
>> 180 muonMinus= G4MuonMinus::MuonMinus();
>> 181
>> 182 mass = charge2 = 0.0;
>> 183
>> 184 nmat = G4Material::GetNumberOfMaterials();
208 const G4MaterialTable* mtable = G4Material:: 185 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
209 if(!pcuts) { pcuts = new G4ProductionCuts(); << 186 pcuts = new G4ProductionCuts();
>> 187
>> 188 couples.resize(nmat,0);
>> 189 cuts.resize(nmat,DBL_MAX);
210 190
211 for(G4int i=0; i<nmat; ++i) { 191 for(G4int i=0; i<nmat; ++i) {
212 couples[i] = new G4MaterialCutsCouple((*mt << 192 couples[i] = new G4MaterialCutsCouple((*mtable)[i],pcuts);
213 } 193 }
214 194
215 dedxElectron = PrepareTable(dedxElectron << 195 splineFlag = G4EmParameters::Instance()->Spline();
216 dedxPositron = PrepareTable(dedxPositron <<
217 dedxMuon = PrepareTable(dedxMuon); <<
218 dedxProton = PrepareTable(dedxProton); <<
219 rangeElectron = PrepareTable(rangeElectro <<
220 rangePositron = PrepareTable(rangePositro <<
221 rangeMuon = PrepareTable(rangeMuon); <<
222 rangeProton = PrepareTable(rangeProton) <<
223 invRangeElectron = PrepareTable(invRangeElec <<
224 invRangePositron = PrepareTable(invRangePosi <<
225 invRangeMuon = PrepareTable(invRangeMuon <<
226 invRangeProton = PrepareTable(invRangeProt <<
227 mscElectron = PrepareTable(mscElectron) <<
228 196
229 builder = new G4LossTableBuilder(true); << 197 dedxElectron = PrepareTable();
230 builder->SetBaseMaterialActive(false); << 198 dedxPositron = PrepareTable();
>> 199 dedxMuon = PrepareTable();
>> 200 dedxProton = PrepareTable();
>> 201 rangeElectron = PrepareTable();
>> 202 rangePositron = PrepareTable();
>> 203 rangeMuon = PrepareTable();
>> 204 rangeProton = PrepareTable();
>> 205 invRangeElectron = PrepareTable();
>> 206 invRangePositron = PrepareTable();
>> 207 invRangeMuon = PrepareTable();
>> 208 invRangeProton = PrepareTable();
>> 209 mscElectron = PrepareTable();
>> 210
>> 211 G4LossTableBuilder builder;
231 212
232 if(verbose>1) { 213 if(verbose>1) {
233 G4cout << "### G4TablesForExtrapolator Bui 214 G4cout << "### G4TablesForExtrapolator Builds electron tables"
234 << G4endl; 215 << G4endl;
235 } 216 }
236 ComputeElectronDEDX(electron, dedxElectron); 217 ComputeElectronDEDX(electron, dedxElectron);
237 builder->BuildRangeTable(dedxElectron,rangeE << 218 builder.BuildRangeTable(dedxElectron,rangeElectron);
238 builder->BuildInverseRangeTable(rangeElectro << 219 builder.BuildInverseRangeTable(rangeElectron, invRangeElectron);
239 220
240 if(verbose>1) { 221 if(verbose>1) {
241 G4cout << "### G4TablesForExtrapolator Bui 222 G4cout << "### G4TablesForExtrapolator Builds positron tables"
242 << G4endl; 223 << G4endl;
243 } 224 }
244 ComputeElectronDEDX(positron, dedxPositron); 225 ComputeElectronDEDX(positron, dedxPositron);
245 builder->BuildRangeTable(dedxPositron, range << 226 builder.BuildRangeTable(dedxPositron, rangePositron);
246 builder->BuildInverseRangeTable(rangePositro << 227 builder.BuildInverseRangeTable(rangePositron, invRangePositron);
247 228
248 if(verbose>1) { 229 if(verbose>1) {
249 G4cout << "### G4TablesForExtrapolator Bui 230 G4cout << "### G4TablesForExtrapolator Builds muon tables" << G4endl;
250 } 231 }
251 ComputeMuonDEDX(muonPlus, dedxMuon); 232 ComputeMuonDEDX(muonPlus, dedxMuon);
252 builder->BuildRangeTable(dedxMuon, rangeMuon << 233 builder.BuildRangeTable(dedxMuon, rangeMuon);
253 builder->BuildInverseRangeTable(rangeMuon, i << 234 builder.BuildInverseRangeTable(rangeMuon, invRangeMuon);
254 << 235 /*
255 if(verbose>2) { << 236 G4cout << "DEDX MUON" << G4endl
256 G4cout << "DEDX MUON" << G4endl; << 237 G4cout << *dedxMuon << G4endl;
257 G4cout << *dedxMuon << G4endl; << 238 G4cout << "RANGE MUON" << G4endl
258 G4cout << "RANGE MUON" << G4endl; << 239 G4cout << *rangeMuon << G4endl;
259 G4cout << *rangeMuon << G4endl; << 240 G4cout << "INVRANGE MUON" << G4endl
260 G4cout << "INVRANGE MUON" << G4endl; << 241 G4cout << *invRangeMuon << G4endl;
261 G4cout << *invRangeMuon << G4endl; << 242 */
262 } <<
263 if(verbose>1) { 243 if(verbose>1) {
264 G4cout << "### G4TablesForExtrapolator Bui 244 G4cout << "### G4TablesForExtrapolator Builds proton tables"
265 << G4endl; 245 << G4endl;
266 } 246 }
267 ComputeProtonDEDX(proton, dedxProton); 247 ComputeProtonDEDX(proton, dedxProton);
268 builder->BuildRangeTable(dedxProton, rangePr << 248 builder.BuildRangeTable(dedxProton, rangeProton);
269 builder->BuildInverseRangeTable(rangeProton, << 249 builder.BuildInverseRangeTable(rangeProton, invRangeProton);
270 250
271 ComputeTrasportXS(electron, mscElectron); 251 ComputeTrasportXS(electron, mscElectron);
272 } 252 }
273 253
274 //....oooOO0OOooo........oooOO0OOooo........oo 254 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
275 255
276 G4PhysicsTable* G4TablesForExtrapolator::Prepa << 256 G4PhysicsTable* G4TablesForExtrapolator::PrepareTable()
277 { 257 {
278 G4PhysicsTable* table = ptr; << 258 G4PhysicsTable* table = new G4PhysicsTable();
279 if(nullptr == ptr) { table = new G4PhysicsTa << 259
280 G4int n = (G4int)table->length(); << 260 for(G4int i=0; i<nmat; ++i) {
281 for(G4int i=n; i<nmat; ++i) { << 261
282 G4PhysicsVector* v = new G4PhysicsLogVecto << 262 G4PhysicsVector* v = new G4PhysicsLogVector(emin, emax, nbins);
>> 263 v->SetSpline(splineFlag);
283 table->push_back(v); 264 table->push_back(v);
284 } 265 }
285 return table; 266 return table;
286 } 267 }
287 268
288 //....oooOO0OOooo........oooOO0OOooo........oo 269 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
289 270
290 void G4TablesForExtrapolator::ComputeElectronD 271 void G4TablesForExtrapolator::ComputeElectronDEDX(
291 const G4Part 272 const G4ParticleDefinition* part,
292 G4PhysicsTable* table) 273 G4PhysicsTable* table)
293 { 274 {
294 G4MollerBhabhaModel* ioni = new G4MollerBhab 275 G4MollerBhabhaModel* ioni = new G4MollerBhabhaModel();
295 G4eBremsstrahlungRelModel* brem = new G4eBre 276 G4eBremsstrahlungRelModel* brem = new G4eBremsstrahlungRelModel();
296 ioni->Initialise(part, cuts); 277 ioni->Initialise(part, cuts);
297 brem->Initialise(part, cuts); 278 brem->Initialise(part, cuts);
298 ioni->SetUseBaseMaterials(false); <<
299 brem->SetUseBaseMaterials(false); <<
300 279
301 mass = electron_mass_c2; 280 mass = electron_mass_c2;
302 charge2 = 1.0; 281 charge2 = 1.0;
303 currentParticle = part; 282 currentParticle = part;
304 283
305 const G4MaterialTable* mtable = G4Material:: 284 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
306 if(0<verbose) { 285 if(0<verbose) {
307 G4cout << "G4TablesForExtrapolator::Comput 286 G4cout << "G4TablesForExtrapolator::ComputeElectronDEDX for "
308 << part->GetParticleName() 287 << part->GetParticleName()
309 << G4endl; 288 << G4endl;
310 } 289 }
311 for(G4int i=0; i<nmat; ++i) { 290 for(G4int i=0; i<nmat; ++i) {
312 291
313 const G4Material* mat = (*mtable)[i]; 292 const G4Material* mat = (*mtable)[i];
314 if(1<verbose) { 293 if(1<verbose) {
315 G4cout << "i= " << i << " mat= " << mat 294 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
316 } 295 }
>> 296 const G4MaterialCutsCouple* couple = couples[i];
317 G4PhysicsVector* aVector = (*table)[i]; 297 G4PhysicsVector* aVector = (*table)[i];
318 298
319 for(G4int j=0; j<=nbins; ++j) { 299 for(G4int j=0; j<=nbins; ++j) {
320 300
321 G4double e = aVector->Energy(j); 301 G4double e = aVector->Energy(j);
322 G4double dedx = ioni->ComputeDEDXPerVol << 302 G4double dedx = ioni->ComputeDEDX(couple,part,e,e)
323 + brem->ComputeDEDXPerVolume(mat,part,e,e); << 303 + brem->ComputeDEDX(couple,part,e,e);
324 if(1<verbose) { 304 if(1<verbose) {
325 G4cout << "j= " << j 305 G4cout << "j= " << j
326 << " e(MeV)= " << e/MeV 306 << " e(MeV)= " << e/MeV
327 << " dedx(Mev/cm)= " << dedx*c 307 << " dedx(Mev/cm)= " << dedx*cm/MeV
328 << " dedx(Mev.cm2/g)= " 308 << " dedx(Mev.cm2/g)= "
329 << dedx/((MeV*mat->GetDensity())/(g/cm2)) 309 << dedx/((MeV*mat->GetDensity())/(g/cm2))
330 << G4endl; 310 << G4endl;
331 } 311 }
332 aVector->PutValue(j,dedx); 312 aVector->PutValue(j,dedx);
333 } 313 }
334 if(splineFlag) { aVector->FillSecondDeriva 314 if(splineFlag) { aVector->FillSecondDerivatives(); }
335 } 315 }
>> 316 delete ioni;
>> 317 delete brem;
336 } 318 }
337 319
338 //....oooOO0OOooo........oooOO0OOooo........oo 320 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
339 321
340 void 322 void
341 G4TablesForExtrapolator::ComputeMuonDEDX(const 323 G4TablesForExtrapolator::ComputeMuonDEDX(const G4ParticleDefinition* part,
342 G4PhysicsTable* table) 324 G4PhysicsTable* table)
343 { 325 {
344 G4BetheBlochModel* ioni = new G4BetheBlochMo 326 G4BetheBlochModel* ioni = new G4BetheBlochModel();
345 G4MuPairProductionModel* pair = new G4MuPair << 327 G4MuPairProductionModel* pair = new G4MuPairProductionModel();
346 G4MuBremsstrahlungModel* brem = new G4MuBrem << 328 G4MuBremsstrahlungModel* brem = new G4MuBremsstrahlungModel();
347 ioni->Initialise(part, cuts); 329 ioni->Initialise(part, cuts);
348 pair->Initialise(part, cuts); 330 pair->Initialise(part, cuts);
349 brem->Initialise(part, cuts); 331 brem->Initialise(part, cuts);
350 ioni->SetUseBaseMaterials(false); <<
351 pair->SetUseBaseMaterials(false); <<
352 brem->SetUseBaseMaterials(false); <<
353 332
354 mass = part->GetPDGMass(); 333 mass = part->GetPDGMass();
355 charge2 = 1.0; 334 charge2 = 1.0;
356 currentParticle = part; 335 currentParticle = part;
357 336
358 const G4MaterialTable* mtable = G4Material:: 337 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
359 338
360 if(0<verbose) { 339 if(0<verbose) {
361 G4cout << "G4TablesForExtrapolator::Comput 340 G4cout << "G4TablesForExtrapolator::ComputeMuonDEDX for "
362 << part->GetParticleName() 341 << part->GetParticleName()
363 << G4endl; 342 << G4endl;
364 } 343 }
365 344
366 for(G4int i=0; i<nmat; ++i) { 345 for(G4int i=0; i<nmat; ++i) {
367 346
368 const G4Material* mat = (*mtable)[i]; 347 const G4Material* mat = (*mtable)[i];
369 if(1<verbose) { 348 if(1<verbose) {
370 G4cout << "i= " << i << " mat= " << mat 349 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
371 } 350 }
>> 351 const G4MaterialCutsCouple* couple = couples[i];
372 G4PhysicsVector* aVector = (*table)[i]; 352 G4PhysicsVector* aVector = (*table)[i];
373 for(G4int j=0; j<=nbins; ++j) { << 353 for(G4int j=0; j<=nbins; j++) {
374 354
375 G4double e = aVector->Energy(j); 355 G4double e = aVector->Energy(j);
376 G4double dedx = ioni->ComputeDEDXPerVol << 356 G4double dedx = ioni->ComputeDEDX(couple,part,e,e) +
377 pair->ComputeDEDXPerVolume(ma << 357 pair->ComputeDEDX(couple,part,e,e) +
378 brem->ComputeDEDXPerVolume(ma << 358 brem->ComputeDEDX(couple,part,e,e);
379 aVector->PutValue(j,dedx); 359 aVector->PutValue(j,dedx);
380 if(1<verbose) { 360 if(1<verbose) {
381 G4cout << "j= " << j 361 G4cout << "j= " << j
382 << " e(MeV)= " << e/MeV 362 << " e(MeV)= " << e/MeV
383 << " dedx(Mev/cm)= " << dedx*c 363 << " dedx(Mev/cm)= " << dedx*cm/MeV
384 << " dedx(Mev/(g/cm2)= " 364 << " dedx(Mev/(g/cm2)= "
385 << dedx/((MeV*mat->GetDensity( 365 << dedx/((MeV*mat->GetDensity())/(g/cm2))
386 << G4endl; 366 << G4endl;
387 } 367 }
388 } 368 }
389 if(splineFlag) { aVector->FillSecondDeriva 369 if(splineFlag) { aVector->FillSecondDerivatives(); }
390 } 370 }
>> 371 delete ioni;
391 } 372 }
392 373
393 //....oooOO0OOooo........oooOO0OOooo........oo 374 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
394 375
395 void 376 void
396 G4TablesForExtrapolator::ComputeProtonDEDX(con 377 G4TablesForExtrapolator::ComputeProtonDEDX(const G4ParticleDefinition* part,
397 G4PhysicsTable* 378 G4PhysicsTable* table)
398 { 379 {
399 G4BetheBlochModel* ioni = new G4BetheBlochMo 380 G4BetheBlochModel* ioni = new G4BetheBlochModel();
400 ioni->Initialise(part, cuts); 381 ioni->Initialise(part, cuts);
401 ioni->SetUseBaseMaterials(false); <<
402 382
403 mass = part->GetPDGMass(); 383 mass = part->GetPDGMass();
404 charge2 = 1.0; 384 charge2 = 1.0;
405 currentParticle = part; 385 currentParticle = part;
406 386
407 const G4MaterialTable* mtable = G4Material:: 387 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
408 388
409 if(0<verbose) { 389 if(0<verbose) {
410 G4cout << "G4TablesForExtrapolator::Comput 390 G4cout << "G4TablesForExtrapolator::ComputeProtonDEDX for "
411 << part->GetParticleName() 391 << part->GetParticleName()
412 << G4endl; 392 << G4endl;
413 } 393 }
414 394
415 for(G4int i=0; i<nmat; ++i) { 395 for(G4int i=0; i<nmat; ++i) {
416 396
417 const G4Material* mat = (*mtable)[i]; 397 const G4Material* mat = (*mtable)[i];
418 if(1<verbose) 398 if(1<verbose)
419 G4cout << "i= " << i << " mat= " << mat 399 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
>> 400 const G4MaterialCutsCouple* couple = couples[i];
420 G4PhysicsVector* aVector = (*table)[i]; 401 G4PhysicsVector* aVector = (*table)[i];
421 for(G4int j=0; j<=nbins; ++j) { << 402 for(G4int j=0; j<=nbins; j++) {
422 403
423 G4double e = aVector->Energy(j); 404 G4double e = aVector->Energy(j);
424 G4double dedx = ioni->ComputeDEDXPerVol << 405 G4double dedx = ioni->ComputeDEDX(couple,part,e,e);
425 aVector->PutValue(j,dedx); 406 aVector->PutValue(j,dedx);
426 if(1<verbose) { 407 if(1<verbose) {
427 G4cout << "j= " << j 408 G4cout << "j= " << j
428 << " e(MeV)= " << e/MeV 409 << " e(MeV)= " << e/MeV
429 << " dedx(Mev/cm)= " << dedx*c 410 << " dedx(Mev/cm)= " << dedx*cm/MeV
430 << " dedx(Mev.cm2/g)= " << ded 411 << " dedx(Mev.cm2/g)= " << dedx/((mat->GetDensity())/(g/cm2))
431 << G4endl; 412 << G4endl;
432 } 413 }
433 } 414 }
434 if(splineFlag) { aVector->FillSecondDeriva 415 if(splineFlag) { aVector->FillSecondDerivatives(); }
435 } 416 }
>> 417 delete ioni;
436 } 418 }
437 419
438 //....oooOO0OOooo........oooOO0OOooo........oo 420 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
439 421
440 void 422 void
441 G4TablesForExtrapolator::ComputeTrasportXS(con 423 G4TablesForExtrapolator::ComputeTrasportXS(const G4ParticleDefinition* part,
442 G4PhysicsTable* table) << 424 G4PhysicsTable* table)
443 { 425 {
444 G4WentzelVIModel* msc = new G4WentzelVIModel 426 G4WentzelVIModel* msc = new G4WentzelVIModel();
445 msc->SetPolarAngleLimit(CLHEP::pi); 427 msc->SetPolarAngleLimit(CLHEP::pi);
446 msc->Initialise(part, cuts); 428 msc->Initialise(part, cuts);
447 msc->SetUseBaseMaterials(false); <<
448 429
449 mass = part->GetPDGMass(); 430 mass = part->GetPDGMass();
450 charge2 = 1.0; 431 charge2 = 1.0;
451 currentParticle = part; 432 currentParticle = part;
452 433
453 const G4MaterialTable* mtable = G4Material:: 434 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
454 435
455 if(0<verbose) { 436 if(0<verbose) {
456 G4cout << "G4TablesForExtrapolator::Comput << 437 G4cout << "G4TablesForExtrapolator::ComputeProtonDEDX for "
457 << part->GetParticleName() 438 << part->GetParticleName()
458 << G4endl; 439 << G4endl;
459 } 440 }
460 441
461 for(G4int i=0; i<nmat; ++i) { 442 for(G4int i=0; i<nmat; ++i) {
462 443
463 const G4Material* mat = (*mtable)[i]; 444 const G4Material* mat = (*mtable)[i];
464 msc->SetCurrentCouple(couples[i]); 445 msc->SetCurrentCouple(couples[i]);
465 if(1<verbose) 446 if(1<verbose)
466 G4cout << "i= " << i << " mat= " << mat 447 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
467 G4PhysicsVector* aVector = (*table)[i]; 448 G4PhysicsVector* aVector = (*table)[i];
468 for(G4int j=0; j<=nbins; ++j) { << 449 for(G4int j=0; j<=nbins; j++) {
469 450
470 G4double e = aVector->Energy(j); 451 G4double e = aVector->Energy(j);
471 G4double xs = msc->CrossSectionPerVolum 452 G4double xs = msc->CrossSectionPerVolume(mat,part,e);
472 aVector->PutValue(j,xs); 453 aVector->PutValue(j,xs);
473 if(1<verbose) { 454 if(1<verbose) {
474 G4cout << "j= " << j << " e(MeV)= " 455 G4cout << "j= " << j << " e(MeV)= " << e/MeV
475 << " xs(1/mm)= " << xs*mm << G 456 << " xs(1/mm)= " << xs*mm << G4endl;
476 } 457 }
477 } 458 }
478 if(splineFlag) { aVector->FillSecondDeriva 459 if(splineFlag) { aVector->FillSecondDerivatives(); }
479 } 460 }
>> 461 delete msc;
480 } 462 }
481 463
482 //....oooOO0OOooo........oooOO0OOooo........oo 464 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
483 465
484 466