| 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: G4LossTableBuilder.cc 92921 2015-09-21 15:06:51Z gcosmo $
>> 27 //
26 // ------------------------------------------- 28 // -------------------------------------------------------------------
27 // 29 //
28 // GEANT4 Class file 30 // GEANT4 Class file
29 // 31 //
30 // 32 //
31 // File name: G4LossTableBuilder 33 // File name: G4LossTableBuilder
32 // 34 //
33 // Author: Vladimir Ivanchenko 35 // Author: Vladimir Ivanchenko
34 // 36 //
35 // Creation date: 03.01.2002 37 // Creation date: 03.01.2002
36 // 38 //
37 // Modifications: 39 // Modifications:
38 // 40 //
39 // 23-01-03 V.Ivanchenko Cut per region 41 // 23-01-03 V.Ivanchenko Cut per region
40 // 21-07-04 V.Ivanchenko Fix problem of range 42 // 21-07-04 V.Ivanchenko Fix problem of range for dedx=0
41 // 08-11-04 Migration to new interface of Stor 43 // 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivanchenko)
42 // 07-12-04 Fix of BuildDEDX table (V.Ivanchen 44 // 07-12-04 Fix of BuildDEDX table (V.Ivanchenko)
43 // 27-03-06 Add bool options isIonisation (V.I 45 // 27-03-06 Add bool options isIonisation (V.Ivanchenko)
44 // 16-01-07 Fill new (not old) DEDX table (V.I 46 // 16-01-07 Fill new (not old) DEDX table (V.Ivanchenko)
45 // 12-02-07 Use G4LPhysicsFreeVector for the i 47 // 12-02-07 Use G4LPhysicsFreeVector for the inverse range table (V.Ivanchenko)
46 // 24-06-09 Removed hidden bin in G4PhysicsVec 48 // 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko)
47 // 49 //
48 // Class Description: 50 // Class Description:
49 // 51 //
50 // ------------------------------------------- 52 // -------------------------------------------------------------------
51 // 53 //
52 //....oooOO0OOooo........oooOO0OOooo........oo 54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
53 //....oooOO0OOooo........oooOO0OOooo........oo 55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
54 56
55 #include "G4LossTableBuilder.hh" 57 #include "G4LossTableBuilder.hh"
56 #include "G4SystemOfUnits.hh" 58 #include "G4SystemOfUnits.hh"
57 #include "G4PhysicsTable.hh" 59 #include "G4PhysicsTable.hh"
58 #include "G4PhysicsLogVector.hh" 60 #include "G4PhysicsLogVector.hh"
59 #include "G4PhysicsTableHelper.hh" 61 #include "G4PhysicsTableHelper.hh"
60 #include "G4PhysicsFreeVector.hh" << 62 #include "G4LPhysicsFreeVector.hh"
61 #include "G4ProductionCutsTable.hh" 63 #include "G4ProductionCutsTable.hh"
62 #include "G4MaterialCutsCouple.hh" 64 #include "G4MaterialCutsCouple.hh"
63 #include "G4Material.hh" 65 #include "G4Material.hh"
64 #include "G4VEmModel.hh" 66 #include "G4VEmModel.hh"
65 #include "G4ParticleDefinition.hh" 67 #include "G4ParticleDefinition.hh"
66 #include "G4LossTableManager.hh" 68 #include "G4LossTableManager.hh"
67 #include "G4EmParameters.hh" <<
68 <<
69 G4bool G4LossTableBuilder::baseMatFlag = false <<
70 std::vector<G4double>* G4LossTableBuilder::the <<
71 std::vector<G4int>* G4LossTableBuilder::the <<
72 std::vector<G4bool>* G4LossTableBuilder::the <<
73 69
74 //....oooOO0OOooo........oooOO0OOooo........oo 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75 71
76 G4LossTableBuilder::G4LossTableBuilder(G4bool << 72 G4LossTableBuilder::G4LossTableBuilder()
77 : isInitializer(master) <<
78 { 73 {
79 theParameters = G4EmParameters::Instance(); << 74 splineFlag = true;
80 if (nullptr == theFlag) { << 75 isInitialized = false;
81 theDensityFactor = new std::vector<G4doubl <<
82 theDensityIdx = new std::vector<G4int>; <<
83 theFlag = new std::vector<G4bool>; <<
84 } <<
85 } <<
86 <<
87 //....oooOO0OOooo........oooOO0OOooo........oo <<
88 76
89 G4LossTableBuilder::~G4LossTableBuilder() << 77 theDensityFactor = new std::vector<G4double>;
90 { << 78 theDensityIdx = new std::vector<G4int>;
91 if (isInitializer) { << 79 theFlag = new std::vector<G4bool>;
92 delete theDensityFactor; <<
93 delete theDensityIdx; <<
94 delete theFlag; <<
95 theDensityFactor = nullptr; <<
96 theDensityIdx = nullptr; <<
97 theFlag = nullptr; <<
98 } <<
99 } <<
100 <<
101 //....oooOO0OOooo........oooOO0OOooo........oo <<
102 <<
103 const std::vector<G4int>* G4LossTableBuilder:: <<
104 { <<
105 return theDensityIdx; <<
106 } <<
107 <<
108 //....oooOO0OOooo........oooOO0OOooo........oo <<
109 <<
110 const std::vector<G4double>* G4LossTableBuilde <<
111 { <<
112 return theDensityFactor; <<
113 } 80 }
114 81
115 //....oooOO0OOooo........oooOO0OOooo........oo 82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
116 83
117 G4bool G4LossTableBuilder::GetFlag(std::size_t << 84 G4LossTableBuilder::~G4LossTableBuilder()
118 { <<
119 return (idx < theFlag->size()) ? (*theFlag)[ <<
120 } <<
121 <<
122 //....oooOO0OOooo........oooOO0OOooo........oo <<
123 <<
124 G4bool G4LossTableBuilder::GetBaseMaterialFlag <<
125 { 85 {
126 return baseMatFlag; << 86 delete theDensityFactor;
>> 87 delete theDensityIdx;
>> 88 delete theFlag;
127 } 89 }
128 90
129 //....oooOO0OOooo........oooOO0OOooo........oo 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
130 92
131 void 93 void
132 G4LossTableBuilder::BuildDEDXTable(G4PhysicsTa 94 G4LossTableBuilder::BuildDEDXTable(G4PhysicsTable* dedxTable,
133 const std:: 95 const std::vector<G4PhysicsTable*>& list)
134 { 96 {
135 InitialiseBaseMaterials(dedxTable); << 97 size_t n_processes = list.size();
136 std::size_t n_processes = list.size(); << 98 //G4cout << "Nproc= " << n_processes << " Ncoup= "
>> 99 //<< dedxTable->size() << G4endl;
137 if(1 >= n_processes) { return; } 100 if(1 >= n_processes) { return; }
138 101
139 std::size_t nCouples = dedxTable->size(); << 102 size_t nCouples = dedxTable->size();
140 //G4cout << "Nproc= " << n_processes << " nC <<
141 // << dedxTable->size() << G4endl; <<
142 if(0 >= nCouples) { return; } 103 if(0 >= nCouples) { return; }
143 104
144 for (std::size_t i=0; i<nCouples; ++i) { << 105 for (size_t i=0; i<nCouples; ++i) {
145 auto pv0 = static_cast<G4PhysicsLogVector* << 106 // if ((*theFlag)[i]) {
146 //if (0 == i) G4cout << i << ". pv0=" << p << 107 G4PhysicsLogVector* pv0 =
147 if(pv0 == nullptr) { continue; } << 108 static_cast<G4PhysicsLogVector*>((*(list[0]))[i]);
148 std::size_t npoints = pv0->GetVectorLength << 109 if(pv0) {
149 auto pv = new G4PhysicsLogVector(*pv0); << 110 size_t npoints = pv0->GetVectorLength();
150 for (std::size_t j=0; j<npoints; ++j) { << 111 G4PhysicsLogVector* pv = new G4PhysicsLogVector(*pv0);
151 G4double dedx = 0.0; << 112 pv->SetSpline(splineFlag);
152 for (std::size_t k=0; k<n_processes; ++k << 113 for (size_t j=0; j<npoints; ++j) {
153 const G4PhysicsVector* pv1 = (*(list[k]))[i] << 114 G4double dedx = 0.0;
154 //if (0 == i) G4cout << " " << k << ". p << 115 for (size_t k=0; k<n_processes; ++k) {
155 dedx += (*pv1)[j]; << 116 G4PhysicsVector* pv1 = (*(list[k]))[i];
>> 117 dedx += (*pv1)[j];
>> 118 }
>> 119 pv->PutValue(j, dedx);
156 } 120 }
157 pv->PutValue(j, dedx); << 121 if(splineFlag) { pv->FillSecondDerivatives(); }
>> 122 G4PhysicsTableHelper::SetPhysicsVector(dedxTable, i, pv);
158 } 123 }
159 if(splineFlag) { pv->FillSecondDerivatives <<
160 G4PhysicsTableHelper::SetPhysicsVector(ded <<
161 } 124 }
162 //G4cout << "### G4LossTableBuilder::BuildDE <<
163 //G4cout << *dedxTable << G4endl; <<
164 } 125 }
165 126
166 //....oooOO0OOooo........oooOO0OOooo........oo 127 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
167 128
168 void G4LossTableBuilder::BuildRangeTable(const 129 void G4LossTableBuilder::BuildRangeTable(const G4PhysicsTable* dedxTable,
169 G4Phy << 130 G4PhysicsTable* rangeTable,
>> 131 G4bool isIonisation)
170 // Build range table from the energy loss tabl 132 // Build range table from the energy loss table
171 { 133 {
172 //G4cout << "### G4LossTableBuilder::BuildRa << 134 size_t nCouples = dedxTable->size();
173 //G4cout << *const_cast<G4PhysicsTable*>(ded <<
174 const std::size_t nCouples = dedxTable->size <<
175 if(0 >= nCouples) { return; } 135 if(0 >= nCouples) { return; }
176 136
177 const std::size_t n = 100; << 137 size_t n = 100;
178 const G4double del = 1.0/(G4double)n; << 138 G4double del = 1.0/(G4double)n;
179 139
180 for (std::size_t i=0; i<nCouples; ++i) { << 140 for (size_t i=0; i<nCouples; ++i) {
181 auto pv = static_cast<G4PhysicsLogVector*> << 141 if(isIonisation) {
182 if((pv == nullptr) || (isBaseMatActive && << 142 if( !(*theFlag)[i] ) { continue; }
183 std::size_t npoints = pv->GetVectorLength( << 143 }
184 std::size_t bin0 = 0; << 144 G4PhysicsLogVector* pv = static_cast<G4PhysicsLogVector*>((*dedxTable)[i]);
>> 145 size_t npoints = pv->GetVectorLength();
>> 146 size_t bin0 = 0;
185 G4double elow = pv->Energy(0); 147 G4double elow = pv->Energy(0);
186 G4double ehigh = pv->Energy(npoints-1); 148 G4double ehigh = pv->Energy(npoints-1);
187 G4double dedx1 = (*pv)[0]; 149 G4double dedx1 = (*pv)[0];
188 150
>> 151 //G4cout << "i= " << i << "npoints= " << npoints << " dedx1= "
>> 152 //<< dedx1 << G4endl;
>> 153
189 // protection for specific cases dedx=0 154 // protection for specific cases dedx=0
190 if(dedx1 == 0.0) { 155 if(dedx1 == 0.0) {
191 for (std::size_t k=1; k<npoints; ++k) { << 156 for (size_t k=1; k<npoints; ++k) {
192 ++bin0; << 157 bin0++;
193 elow = pv->Energy(k); << 158 elow = pv->Energy(k);
194 dedx1 = (*pv)[k]; << 159 dedx1 = (*pv)[k];
195 if(dedx1 > 0.0) { break; } << 160 if(dedx1 > 0.0) { break; }
196 } 161 }
197 npoints -= bin0; 162 npoints -= bin0;
198 } 163 }
>> 164 //G4cout<<"New Range vector" << G4endl;
>> 165 //G4cout<<"nbins= "<<npoints-1<<" elow= "<<elow<<" ehigh= "<<ehigh
>> 166 // <<" bin0= " << bin0 <<G4endl;
199 167
200 // initialisation of a new vector 168 // initialisation of a new vector
201 if(npoints < 3) { npoints = 3; } << 169 if(npoints < 2) { npoints = 2; }
202 170
203 delete (*rangeTable)[i]; 171 delete (*rangeTable)[i];
204 G4PhysicsLogVector* v; 172 G4PhysicsLogVector* v;
205 if(0 == bin0) { v = new G4PhysicsLogVector 173 if(0 == bin0) { v = new G4PhysicsLogVector(*pv); }
206 else { v = new G4PhysicsLogVector(elow, eh << 174 else { v = new G4PhysicsLogVector(elow, ehigh, npoints-1); }
>> 175
>> 176 // dedx is exact zero cannot build range table
>> 177 if(2 == npoints) {
>> 178 v->PutValue(0,1000.);
>> 179 v->PutValue(1,2000.);
>> 180 G4PhysicsTableHelper::SetPhysicsVector(rangeTable, i, v);
>> 181 return;
>> 182 }
>> 183 v->SetSpline(splineFlag);
207 184
208 // assumed dedx proportional to beta 185 // assumed dedx proportional to beta
209 G4double energy1 = v->Energy(0); 186 G4double energy1 = v->Energy(0);
210 G4double range = 2.*energy1/dedx1; 187 G4double range = 2.*energy1/dedx1;
211 /* << 188 //G4cout << "range0= " << range << G4endl;
212 G4cout << "New Range vector Npoints=" << v <<
213 << " coupleIdx=" << i << " spline=" << v- <<
214 << " Elow=" << v->GetMinEnergy() <<" Ehig <<
215 << " DEDX(Elow)=" << dedx1 << " R(Elow)=" <<
216 */ <<
217 v->PutValue(0,range); 189 v->PutValue(0,range);
218 190
219 for (std::size_t j=1; j<npoints; ++j) { << 191 for (size_t j=1; j<npoints; ++j) {
220 192
221 G4double energy2 = v->Energy(j); 193 G4double energy2 = v->Energy(j);
222 G4double de = (energy2 - energy1) * 194 G4double de = (energy2 - energy1) * del;
223 G4double energy = energy2 + de*0.5; 195 G4double energy = energy2 + de*0.5;
224 G4double sum = 0.0; 196 G4double sum = 0.0;
225 std::size_t idx = j - 1; << 197 //G4cout << "j= " << j << " e1= " << energy1 << " e2= " << energy2
226 for (std::size_t k=0; k<n; ++k) { << 198 // << " n= " << n << G4endl;
227 energy -= de; << 199 for (size_t k=0; k<n; ++k) {
228 dedx1 = pv->Value(energy, idx); << 200 energy -= de;
229 if(dedx1 > 0.0) { sum += de/dedx1; } << 201 dedx1 = pv->Value(energy);
>> 202 if(dedx1 > 0.0) { sum += de/dedx1; }
230 } 203 }
231 range += sum; 204 range += sum;
232 /* <<
233 if(energy < 10.) <<
234 G4cout << "j= " << j << " e1= " << energy1 < <<
235 << " n= " << n << " range=" << range< <<
236 */ <<
237 v->PutValue(j,range); 205 v->PutValue(j,range);
238 energy1 = energy2; 206 energy1 = energy2;
239 } 207 }
240 if(splineFlag) { v->FillSecondDerivatives( 208 if(splineFlag) { v->FillSecondDerivatives(); }
241 G4PhysicsTableHelper::SetPhysicsVector(ran 209 G4PhysicsTableHelper::SetPhysicsVector(rangeTable, i, v);
242 } 210 }
243 //G4cout << "### Range table" << G4endl; <<
244 //G4cout << *rangeTable << G4endl; <<
245 } 211 }
246 212
247 //....oooOO0OOooo........oooOO0OOooo........oo 213 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
248 214
249 void 215 void
250 G4LossTableBuilder::BuildInverseRangeTable(con 216 G4LossTableBuilder::BuildInverseRangeTable(const G4PhysicsTable* rangeTable,
251 G4P << 217 G4PhysicsTable* invRangeTable,
>> 218 G4bool isIonisation)
252 // Build inverse range table from the energy l 219 // Build inverse range table from the energy loss table
253 { 220 {
254 std::size_t nCouples = rangeTable->size(); << 221 size_t nCouples = rangeTable->size();
255 if(0 >= nCouples) { return; } 222 if(0 >= nCouples) { return; }
256 223
257 for (std::size_t i=0; i<nCouples; ++i) { << 224 for (size_t i=0; i<nCouples; ++i) {
>> 225
>> 226 if(isIonisation) {
>> 227 if( !(*theFlag)[i] ) { continue; }
>> 228 }
258 G4PhysicsVector* pv = (*rangeTable)[i]; 229 G4PhysicsVector* pv = (*rangeTable)[i];
259 if((pv == nullptr) || (isBaseMatActive && << 230 size_t npoints = pv->GetVectorLength();
260 std::size_t npoints = pv->GetVectorLength( << 231 G4double rlow = (*pv)[0];
>> 232 G4double rhigh = (*pv)[npoints-1];
261 233
262 delete (*invRangeTable)[i]; 234 delete (*invRangeTable)[i];
263 auto v = new G4PhysicsFreeVector(npoints, << 235 G4LPhysicsFreeVector* v = new G4LPhysicsFreeVector(npoints,rlow,rhigh);
>> 236 v->SetSpline(splineFlag);
264 237
265 for (std::size_t j=0; j<npoints; ++j) { << 238 for (size_t j=0; j<npoints; ++j) {
266 G4double e = pv->Energy(j); 239 G4double e = pv->Energy(j);
267 G4double r = (*pv)[j]; 240 G4double r = (*pv)[j];
268 v->PutValues(j,r,e); 241 v->PutValues(j,r,e);
269 } 242 }
270 if (splineFlag) { v->FillSecondDerivatives << 243 if(splineFlag) { v->FillSecondDerivatives(); }
271 v->EnableLogBinSearch(theParameters->Numbe <<
272 244
273 G4PhysicsTableHelper::SetPhysicsVector(inv 245 G4PhysicsTableHelper::SetPhysicsVector(invRangeTable, i, v);
274 } 246 }
275 //G4cout << "### Inverse range table" << G4e <<
276 //G4cout << *invRangeTable << G4endl; <<
277 } 247 }
278 248
279 //....oooOO0OOooo........oooOO0OOooo........oo 249 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
280 250
281 void G4LossTableBuilder::InitialiseBaseMateria << 251 void
>> 252 G4LossTableBuilder::InitialiseBaseMaterials(G4PhysicsTable* table)
282 { 253 {
283 if(!isInitializer) { return; } << 254 size_t nCouples = table->size();
>> 255 size_t nFlags = theFlag->size();
>> 256
>> 257 if(nCouples == nFlags && isInitialized) { return; }
>> 258
>> 259 isInitialized = true;
>> 260
>> 261 //G4cout << "%%%%%% G4LossTableBuilder::InitialiseBaseMaterials Ncouples= "
>> 262 // << nCouples << " FlagSize= " << nFlags << G4endl;
>> 263
>> 264 // variable density check
284 const G4ProductionCutsTable* theCoupleTable= 265 const G4ProductionCutsTable* theCoupleTable=
285 G4ProductionCutsTable::GetProductionCutsTa 266 G4ProductionCutsTable::GetProductionCutsTable();
286 std::size_t nCouples = theCoupleTable->GetTa << 267
287 std::size_t nFlags = theFlag->size(); <<
288 /* 268 /*
289 G4cout << "### InitialiseBaseMaterials: nCou << 269 for(size_t i=0; i<nFlags; ++i) {
290 << " nFlags=" << nFlags << " isInit:" << is << 270 G4cout << "CoupleIdx= " << i << " Flag= " << (*theFlag)[i]
291 << " baseMat:" << baseMatFlag << G4endl; << 271 << " tableFlag= " << table->GetFlag(i) << " "
>> 272 << theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial()->GetName()
>> 273 << G4endl;
>> 274 }
292 */ 275 */
293 // define base material flag << 276
294 if(isBaseMatActive && !baseMatFlag) { << 277 // expand vectors
295 for(G4int i=0; i<(G4int)nCouples; ++i) { << 278 if(nFlags < nCouples) {
296 if(nullptr != theCoupleTable->GetMateria << 279 for(size_t i=nFlags; i<nCouples; ++i) {
297 baseMatFlag = true; << 280 theDensityFactor->push_back(1.0);
298 isInitialized = false; << 281 }
299 break; << 282 for(size_t i=nFlags; i<nCouples; ++i) { theDensityIdx->push_back(-1); }
>> 283 for(size_t i=nFlags; i<nCouples; ++i) { theFlag->push_back(true); }
>> 284 }
>> 285 for(size_t i=0; i<nCouples; ++i) {
>> 286
>> 287 // base material is needed only for a couple which is not
>> 288 // initialised and for which tables will be computed
>> 289 (*theFlag)[i] = table->GetFlag(i);
>> 290 if ((*theDensityIdx)[i] < 0) {
>> 291 (*theDensityIdx)[i] = i;
>> 292 const G4MaterialCutsCouple* couple =
>> 293 theCoupleTable->GetMaterialCutsCouple(i);
>> 294 const G4ProductionCuts* pcuts = couple->GetProductionCuts();
>> 295 const G4Material* mat = couple->GetMaterial();
>> 296 const G4Material* bmat = mat->GetBaseMaterial();
>> 297
>> 298 // base material exists - find it and check if it can be reused
>> 299 if(bmat) {
>> 300 for(size_t j=0; j<nCouples; ++j) {
>> 301
>> 302 if(j == i) { continue; }
>> 303 const G4MaterialCutsCouple* bcouple =
>> 304 theCoupleTable->GetMaterialCutsCouple(j);
>> 305
>> 306 if(bcouple->GetMaterial() == bmat &&
>> 307 bcouple->GetProductionCuts() == pcuts) {
>> 308
>> 309 // based couple exist in the same region
>> 310 (*theDensityIdx)[i] = j;
>> 311 (*theDensityFactor)[i] = mat->GetDensity()/bmat->GetDensity();
>> 312 (*theFlag)[i] = false;
>> 313
>> 314 // ensure that there will no double initialisation
>> 315 (*theDensityIdx)[j] = j;
>> 316 (*theDensityFactor)[j] = 1.0;
>> 317 (*theFlag)[j] = true;
>> 318 break;
>> 319 }
>> 320 }
300 } 321 }
301 } 322 }
302 } 323 }
>> 324 /*
>> 325 for(size_t i=0; i<nCouples; ++i) {
>> 326 G4cout << "CoupleIdx= " << i << " Flag= " << (*theFlag)[i]
>> 327 << " TableFlag= " << table->GetFlag(i) << " "
>> 328 << theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial()->GetName()
>> 329 << G4endl;
>> 330 }
>> 331 G4cout << "%%%%%% G4LossTableBuilder::InitialiseBaseMaterials end"
>> 332 << G4endl;
>> 333 */
>> 334 }
>> 335
>> 336 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 337
>> 338 void G4LossTableBuilder::InitialiseCouples()
>> 339 {
>> 340 isInitialized = true;
303 341
304 if(nFlags != nCouples) { isInitialized = fal << 342 // variable density initialisation for the cas without tables
305 if(isInitialized) { return; } << 343 const G4ProductionCutsTable* theCoupleTable=
>> 344 G4ProductionCutsTable::GetProductionCutsTable();
>> 345 size_t nCouples = theCoupleTable->GetTableSize();
>> 346 //G4cout << "%%%%%% G4LossTableBuilder::InitialiseCouples() nCouples= "
>> 347 // << nCouples << G4endl;
306 348
307 // reserve memory << 349 theDensityFactor->resize(nCouples, 1.0);
>> 350 theDensityIdx->resize(nCouples, -1);
308 theFlag->resize(nCouples, true); 351 theFlag->resize(nCouples, true);
309 theDensityFactor->resize(nCouples,1.0); <<
310 theDensityIdx->resize(nCouples, 0); <<
311 352
312 // define default flag and index of used mat << 353 for(size_t i=0; i<nCouples; ++i) {
313 for (G4int i=0; i<(G4int)nCouples; ++i) { <<
314 (*theFlag)[i] = (nullptr == table) ? true <<
315 (*theDensityIdx)[i] = i; <<
316 } <<
317 isInitialized = true; <<
318 if (!baseMatFlag) { return; } <<
319 354
320 // use base materials <<
321 for (G4int i=0; i<(G4int)nCouples; ++i) { <<
322 // base material is needed only for a coup 355 // base material is needed only for a couple which is not
323 // initialised and for which tables will b 356 // initialised and for which tables will be computed
324 auto couple = theCoupleTable->GetMaterialC << 357 if ((*theDensityIdx)[i] < 0) {
325 auto pcuts = couple->GetProductionCuts(); << 358 (*theDensityIdx)[i] = i;
326 auto mat = couple->GetMaterial(); << 359 const G4MaterialCutsCouple* couple =
327 auto bmat = mat->GetBaseMaterial(); << 360 theCoupleTable->GetMaterialCutsCouple(i);
328 << 361 const G4ProductionCuts* pcuts = couple->GetProductionCuts();
329 // base material exists - find it and chec << 362 const G4Material* mat = couple->GetMaterial();
330 if(nullptr != bmat) { << 363 const G4Material* bmat = mat->GetBaseMaterial();
331 for(G4int j=0; j<(G4int)nCouples; ++j) { << 364
332 if(j == i) { continue; } << 365 // base material exists - find it and check if it can be reused
333 auto bcouple = theCoupleTable->GetMaterialCu << 366 if(bmat) {
334 << 367 for(size_t j=0; j<nCouples; ++j) {
335 if(bcouple->GetMaterial() == bmat && << 368
336 bcouple->GetProductionCuts() == pcuts) { << 369 if(j == i) { continue; }
337 << 370 const G4MaterialCutsCouple* bcouple =
338 // based couple exist in the same region << 371 theCoupleTable->GetMaterialCutsCouple(j);
339 (*theDensityFactor)[i] = mat->GetDensity() << 372
340 (*theDensityIdx)[i] = j; << 373 if(bcouple->GetMaterial() == bmat &&
341 (*theFlag)[i] = false; << 374 bcouple->GetProductionCuts() == pcuts) {
342 << 375
343 // ensure that there will no double initia << 376 // based couple exist in the same region
344 (*theDensityFactor)[j] = 1.0; << 377 (*theDensityIdx)[i] = j;
345 (*theDensityIdx)[j] = j; << 378 (*theDensityFactor)[i] = mat->GetDensity()/bmat->GetDensity();
346 (*theFlag)[j] = true; << 379 (*theFlag)[i] = false;
347 break; << 380
348 } << 381 // ensure that there will no double initialisation
>> 382 (*theDensityIdx)[j] = j;
>> 383 (*theDensityFactor)[j] = 1.0;
>> 384 (*theFlag)[j] = true;
>> 385 break;
>> 386 }
>> 387 }
349 } 388 }
350 } 389 }
351 } 390 }
>> 391 /*
>> 392 for(size_t i=0; i<nCouples; ++i) {
>> 393 G4cout << "CoupleIdx= " << i << " Flag= " << (*theFlag)[i] << " "
>> 394 << theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial()->GetName()
>> 395 << " DensityFactor= " << (*theDensityFactor)[i]
>> 396 << G4endl;
>> 397 }
>> 398 G4cout << "%%%%%% G4LossTableBuilder::InitialiseCouples() end" << G4endl;
>> 399 */
352 } 400 }
353 401
354 //....oooOO0OOooo........oooOO0OOooo........oo 402 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
355 403
356 G4PhysicsTable* 404 G4PhysicsTable*
357 G4LossTableBuilder::BuildTableForModel(G4Physi 405 G4LossTableBuilder::BuildTableForModel(G4PhysicsTable* aTable,
358 G4VEmMo 406 G4VEmModel* model,
359 const G 407 const G4ParticleDefinition* part,
360 G4doubl 408 G4double emin, G4double emax,
361 G4bool 409 G4bool spline)
362 { 410 {
363 // check input 411 // check input
364 G4PhysicsTable* table = G4PhysicsTableHelper 412 G4PhysicsTable* table = G4PhysicsTableHelper::PreparePhysicsTable(aTable);
365 if (nullptr == table) { return table; } << 413 if(!table) { return table; }
366 if (aTable != nullptr && aTable != table) { << 414 if(emin >= emax) {
367 aTable->clearAndDestroy(); << 415 table->clearAndDestroy();
368 delete aTable; << 416 delete table;
>> 417 table = 0;
>> 418 return table;
369 } 419 }
370 <<
371 InitialiseBaseMaterials(table); 420 InitialiseBaseMaterials(table);
372 G4int nbins = theParameters->NumberOfBinsPer << 421
>> 422 G4int nbins = G4LossTableManager::Instance()->GetNumberOfBinsPerDecade();
373 423
374 // Access to materials 424 // Access to materials
375 const G4ProductionCutsTable* theCoupleTable= 425 const G4ProductionCutsTable* theCoupleTable=
376 G4ProductionCutsTable::GetProductionCu 426 G4ProductionCutsTable::GetProductionCutsTable();
377 std::size_t numOfCouples = theCoupleTable->G << 427 size_t numOfCouples = theCoupleTable->GetTableSize();
378 /* << 428
379 G4cout << " G4LossTableBuilder::BuildTable << 429 G4PhysicsLogVector* aVector = 0;
380 << " isMaster=" << isInitializer << " model <<
381 << " " << part->GetParticleName() << G4end <<
382 */ <<
383 G4PhysicsLogVector* aVector = nullptr; <<
384 430
385 for(G4int i=0; i<(G4int)numOfCouples; ++i) { << 431 for(size_t i=0; i<numOfCouples; ++i) {
386 //G4cout << i << ". " << (*theFlag)[i] << << 432
387 if (table->GetFlag(i)) { << 433 //G4cout<< "i= " << i << " Flag= " << GetFlag(i) << G4endl;
>> 434
>> 435 if (GetFlag(i)) {
388 436
389 // create physics vector and fill it 437 // create physics vector and fill it
390 auto couple = theCoupleTable->GetMateria << 438 const G4MaterialCutsCouple* couple =
>> 439 theCoupleTable->GetMaterialCutsCouple(i);
391 delete (*table)[i]; 440 delete (*table)[i];
392 441
393 // if start from zero then change the sc 442 // if start from zero then change the scale
>> 443
394 const G4Material* mat = couple->GetMater 444 const G4Material* mat = couple->GetMaterial();
395 445
396 G4double tmin = std::max(emin, model->Mi << 446 G4double tmin = std::max(emin,model->MinPrimaryEnergy(mat,part));
397 if(0.0 >= tmin) { tmin = CLHEP::eV; } << 447 if(0.0 >= tmin) { tmin = eV; }
398 G4int n = nbins; 448 G4int n = nbins;
399 449
400 if(tmin >= emax) { 450 if(tmin >= emax) {
401 aVector = nullptr; << 451 aVector = 0;
402 } else { 452 } else {
403 n *= G4lrint(std::log10(emax/tmin)); << 453 n *= G4int(std::log10(emax/tmin) + 0.5);
404 n = std::max(n, 3); << 454 if(n < 3) { n = 3; }
405 aVector = new G4PhysicsLogVector(tmin, << 455 aVector = new G4PhysicsLogVector(tmin, emax, n);
406 } 456 }
407 457
408 if(nullptr != aVector) { << 458 if(aVector) {
409 //G4cout << part->GetParticleName() << << 459 aVector->SetSpline(spline);
410 // << " emin= " << tmin << " emax=" << e <<
411 for(G4int j=0; j<=n; ++j) { 460 for(G4int j=0; j<=n; ++j) {
412 G4double e = aVector->Energy(j); << 461 aVector->PutValue(j, model->Value(couple, part,
413 G4double y = model->Value(couple, part, e) << 462 aVector->Energy(j)));
414 //G4cout << " " << j << ") E=" << e < <<
415 aVector->PutValue(j, y); <<
416 } 463 }
417 if(spline) { aVector->FillSecondDeriva 464 if(spline) { aVector->FillSecondDerivatives(); }
418 } 465 }
419 G4PhysicsTableHelper::SetPhysicsVector(t 466 G4PhysicsTableHelper::SetPhysicsVector(table, i, aVector);
420 } 467 }
421 } 468 }
422 /* 469 /*
423 G4cout << "G4LossTableBuilder::BuildTableFor 470 G4cout << "G4LossTableBuilder::BuildTableForModel done for "
424 << part->GetParticleName() << " and " 471 << part->GetParticleName() << " and "<< model->GetName()
425 << " " << table << G4endl; 472 << " " << table << G4endl;
426 */ 473 */
427 //G4cout << *table << G4endl; <<
428 return table; 474 return table;
429 } 475 }
430 476
431 //....oooOO0OOooo........oooOO0OOooo........oo 477 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 478
432 479