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