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