| 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 // G4ProductionCutsTable class implementation 26 // G4ProductionCutsTable class implementation
27 // 27 //
28 // Author: M.Asai, 5 October 2002 - First impl 28 // Author: M.Asai, 5 October 2002 - First implementation
29 // Modifications: H.Kurashige, 2004-2008 29 // Modifications: H.Kurashige, 2004-2008
30 // ------------------------------------------- 30 // --------------------------------------------------------------------
31 31
32 #include "G4ProductionCutsTable.hh" 32 #include "G4ProductionCutsTable.hh"
33 #include "G4ProductionCuts.hh" 33 #include "G4ProductionCuts.hh"
34 #include "G4MCCIndexConversionTable.hh" 34 #include "G4MCCIndexConversionTable.hh"
35 #include "G4ProductionCutsTableMessenger.hh" 35 #include "G4ProductionCutsTableMessenger.hh"
36 #include "G4ParticleDefinition.hh" 36 #include "G4ParticleDefinition.hh"
37 #include "G4ParticleTable.hh" 37 #include "G4ParticleTable.hh"
38 #include "G4RegionStore.hh" 38 #include "G4RegionStore.hh"
39 #include "G4LogicalVolume.hh" 39 #include "G4LogicalVolume.hh"
40 #include "G4VPhysicalVolume.hh" 40 #include "G4VPhysicalVolume.hh"
41 #include "G4RToEConvForElectron.hh" 41 #include "G4RToEConvForElectron.hh"
42 #include "G4RToEConvForGamma.hh" 42 #include "G4RToEConvForGamma.hh"
43 #include "G4RToEConvForPositron.hh" 43 #include "G4RToEConvForPositron.hh"
44 #include "G4RToEConvForProton.hh" 44 #include "G4RToEConvForProton.hh"
45 #include "G4MaterialTable.hh" 45 #include "G4MaterialTable.hh"
46 #include "G4Material.hh" 46 #include "G4Material.hh"
47 #include "G4UnitsTable.hh" 47 #include "G4UnitsTable.hh"
48 48
49 #include "G4Timer.hh" 49 #include "G4Timer.hh"
50 #include "G4SystemOfUnits.hh" 50 #include "G4SystemOfUnits.hh"
51 #include "G4ios.hh" 51 #include "G4ios.hh"
52 #include <iomanip> 52 #include <iomanip>
53 #include <fstream> 53 #include <fstream>
54 54
55 G4ProductionCutsTable* G4ProductionCutsTable:: 55 G4ProductionCutsTable* G4ProductionCutsTable::fProductionCutsTable = nullptr;
56 56
57 // ------------------------------------------- 57 // --------------------------------------------------------------------
58 G4ProductionCutsTable* G4ProductionCutsTable:: 58 G4ProductionCutsTable* G4ProductionCutsTable::GetProductionCutsTable()
59 { 59 {
60 static G4ProductionCutsTable theProductionC 60 static G4ProductionCutsTable theProductionCutsTable;
61 if(fProductionCutsTable == nullptr) 61 if(fProductionCutsTable == nullptr)
62 { 62 {
63 fProductionCutsTable = &theProductionCuts 63 fProductionCutsTable = &theProductionCutsTable;
64 } 64 }
65 return fProductionCutsTable; 65 return fProductionCutsTable;
66 } 66 }
67 67
68 // ------------------------------------------- 68 // --------------------------------------------------------------------
69 G4ProductionCutsTable::G4ProductionCutsTable() 69 G4ProductionCutsTable::G4ProductionCutsTable()
70 { 70 {
71 for(std::size_t i=0; i< NumberOfG4CutIndex; 71 for(std::size_t i=0; i< NumberOfG4CutIndex; ++i)
72 { 72 {
73 rangeCutTable.push_back(new std::vector<G4 73 rangeCutTable.push_back(new std::vector<G4double>);
74 energyCutTable.push_back(new std::vector<G 74 energyCutTable.push_back(new std::vector<G4double>);
75 rangeDoubleVector[i] = nullptr; 75 rangeDoubleVector[i] = nullptr;
76 energyDoubleVector[i] = nullptr; 76 energyDoubleVector[i] = nullptr;
77 converters[i] = nullptr; 77 converters[i] = nullptr;
78 } 78 }
79 fG4RegionStore = G4RegionStore::GetInstance( 79 fG4RegionStore = G4RegionStore::GetInstance();
80 defaultProductionCuts = new G4ProductionCuts 80 defaultProductionCuts = new G4ProductionCuts();
81 81
82 // add messenger for UI 82 // add messenger for UI
83 fMessenger = new G4ProductionCutsTableMessen 83 fMessenger = new G4ProductionCutsTableMessenger(this);
84 } 84 }
85 85
86 // ------------------------------------------- 86 // --------------------------------------------------------------------
87 G4ProductionCutsTable::~G4ProductionCutsTable( 87 G4ProductionCutsTable::~G4ProductionCutsTable()
88 { 88 {
89 delete defaultProductionCuts; 89 delete defaultProductionCuts;
90 defaultProductionCuts = nullptr; 90 defaultProductionCuts = nullptr;
91 91
92 for(auto itr=coupleTable.cbegin(); itr!=coup 92 for(auto itr=coupleTable.cbegin(); itr!=coupleTable.cend(); ++itr)
93 { 93 {
94 delete (*itr); 94 delete (*itr);
95 } 95 }
96 coupleTable.clear(); 96 coupleTable.clear();
97 97
98 for(std::size_t i=0; i< NumberOfG4CutIndex; 98 for(std::size_t i=0; i< NumberOfG4CutIndex; ++i)
99 { 99 {
100 delete rangeCutTable[i]; 100 delete rangeCutTable[i];
101 delete energyCutTable[i]; 101 delete energyCutTable[i];
102 delete converters[i]; 102 delete converters[i];
103 if(rangeDoubleVector[i] != nullptr) delete 103 if(rangeDoubleVector[i] != nullptr) delete [] rangeDoubleVector[i];
104 if(energyDoubleVector[i] != nullptr) delet 104 if(energyDoubleVector[i] != nullptr) delete [] energyDoubleVector[i];
105 rangeCutTable[i] = nullptr; 105 rangeCutTable[i] = nullptr;
106 energyCutTable[i] = nullptr; 106 energyCutTable[i] = nullptr;
107 converters[i] = nullptr; 107 converters[i] = nullptr;
108 rangeDoubleVector[i] = nullptr; 108 rangeDoubleVector[i] = nullptr;
109 energyDoubleVector[i] = nullptr; 109 energyDoubleVector[i] = nullptr;
110 if(i < 4) <<
111 { <<
112 delete userEnergyCuts[i]; <<
113 } <<
114 } 110 }
115 fProductionCutsTable = nullptr; 111 fProductionCutsTable = nullptr;
116 112
117 delete fMessenger; 113 delete fMessenger;
118 fMessenger = nullptr; 114 fMessenger = nullptr;
119 } 115 }
120 116
121 // ------------------------------------------- 117 // --------------------------------------------------------------------
122 void G4ProductionCutsTable::SetEnergyCutVector << 118 void G4ProductionCutsTable::UpdateCoupleTable(G4VPhysicalVolume* /*currWorld*/)
123 <<
124 { 119 {
125 if (idx >= 4) { << 120 if(firstUse)
126 G4ExceptionDescription ed; << 121 {
127 ed << "Wrong index= " << idx << "; it shou << 122 if(G4ParticleTable::GetParticleTable()->FindParticle("gamma"))
128 G4Exception("G4ProductionCutsTable::SetEne << 123 {
129 "CUTS0100", FatalException, ed); << 124 converters[0] = new G4RToEConvForGamma();
130 return; << 125 converters[0]->SetVerboseLevel(GetVerboseLevel());
>> 126 }
>> 127 if(G4ParticleTable::GetParticleTable()->FindParticle("e-"))
>> 128 {
>> 129 converters[1] = new G4RToEConvForElectron();
>> 130 converters[1]->SetVerboseLevel(GetVerboseLevel());
>> 131 }
>> 132 if(G4ParticleTable::GetParticleTable()->FindParticle("e+"))
>> 133 {
>> 134 converters[2] = new G4RToEConvForPositron();
>> 135 converters[2]->SetVerboseLevel(GetVerboseLevel());
>> 136 }
>> 137 if(G4ParticleTable::GetParticleTable()->FindParticle("proton"))
>> 138 {
>> 139 converters[3] = new G4RToEConvForProton();
>> 140 converters[3]->SetVerboseLevel(GetVerboseLevel());
>> 141 }
>> 142 firstUse = false;
131 } 143 }
132 userEnergyCuts[idx] = new std::vector<G4doub <<
133 } <<
134 144
135 // ------------------------------------------- <<
136 void G4ProductionCutsTable::CreateCoupleTables <<
137 { <<
138 // Reset "used" flags of all couples 145 // Reset "used" flags of all couples
139 for(auto CoupleItr=coupleTable.cbegin(); 146 for(auto CoupleItr=coupleTable.cbegin();
140 CoupleItr!=coupleTable.cend(); ++Co 147 CoupleItr!=coupleTable.cend(); ++CoupleItr)
141 { 148 {
142 (*CoupleItr)->SetUseFlag(false); 149 (*CoupleItr)->SetUseFlag(false);
143 } 150 }
144 151
145 // Update Material-Cut-Couple 152 // Update Material-Cut-Couple
146 for(auto rItr=fG4RegionStore->cbegin(); rItr 153 for(auto rItr=fG4RegionStore->cbegin(); rItr!=fG4RegionStore->cend(); ++rItr)
147 { 154 {
148 // Material scan is to be done only for th 155 // Material scan is to be done only for the regions appear in the
149 // current tracking world. 156 // current tracking world.
150 // if((*rItr)->GetWorldPhysical()!=curr 157 // if((*rItr)->GetWorldPhysical()!=currentWorld) continue;
151 158
152 if( (*rItr)->IsInMassGeometry() || (*rItr) 159 if( (*rItr)->IsInMassGeometry() || (*rItr)->IsInParallelGeometry() )
153 { 160 {
154 G4ProductionCuts* fProductionCut = (*rIt 161 G4ProductionCuts* fProductionCut = (*rItr)->GetProductionCuts();
155 auto mItr = (*rItr)->GetMaterialIterator 162 auto mItr = (*rItr)->GetMaterialIterator();
156 std::size_t nMaterial = (*rItr)->GetNumb 163 std::size_t nMaterial = (*rItr)->GetNumberOfMaterials();
157 (*rItr)->ClearMap(); 164 (*rItr)->ClearMap();
158 165
159 for(std::size_t iMate=0; iMate<nMaterial 166 for(std::size_t iMate=0; iMate<nMaterial; ++iMate)
160 { 167 {
161 //check if this material cut couple ha 168 //check if this material cut couple has already been made
162 G4bool coupleAlreadyDefined = false; 169 G4bool coupleAlreadyDefined = false;
163 G4MaterialCutsCouple* aCouple; 170 G4MaterialCutsCouple* aCouple;
164 for(auto cItr=coupleTable.cbegin(); cI 171 for(auto cItr=coupleTable.cbegin(); cItr!=coupleTable.cend(); ++cItr)
165 { 172 {
166 if( (*cItr)->GetMaterial()==(*mItr) 173 if( (*cItr)->GetMaterial()==(*mItr)
167 && (*cItr)->GetProductionCuts()==fP 174 && (*cItr)->GetProductionCuts()==fProductionCut)
168 { 175 {
169 coupleAlreadyDefined = true; 176 coupleAlreadyDefined = true;
170 aCouple = *cItr; 177 aCouple = *cItr;
171 break; 178 break;
172 } 179 }
173 } 180 }
174 181
175 // If this combination is new, cleate 182 // If this combination is new, cleate and register a couple
176 if(!coupleAlreadyDefined) 183 if(!coupleAlreadyDefined)
177 { 184 {
178 aCouple = new G4MaterialCutsCouple(( 185 aCouple = new G4MaterialCutsCouple((*mItr),fProductionCut);
179 coupleTable.push_back(aCouple); 186 coupleTable.push_back(aCouple);
180 aCouple->SetIndex(G4int(coupleTable. << 187 aCouple->SetIndex(coupleTable.size()-1);
181 } 188 }
182 189
183 // Register this couple to the region 190 // Register this couple to the region
184 (*rItr)->RegisterMaterialCouplePair((* 191 (*rItr)->RegisterMaterialCouplePair((*mItr),aCouple);
185 192
186 // Set the couple to the proper logica 193 // Set the couple to the proper logical volumes in that region
187 aCouple->SetUseFlag(); 194 aCouple->SetUseFlag();
188 195
189 auto rootLVItr = (*rItr)->GetRootLogic 196 auto rootLVItr = (*rItr)->GetRootLogicalVolumeIterator();
190 std::size_t nRootLV = (*rItr)->GetNumb 197 std::size_t nRootLV = (*rItr)->GetNumberOfRootVolumes();
191 for(std::size_t iLV=0; iLV<nRootLV; ++ 198 for(std::size_t iLV=0; iLV<nRootLV; ++iLV)
192 { 199 {
193 // Set the couple to the proper logi 200 // Set the couple to the proper logical volumes in that region
194 G4LogicalVolume* aLV = *rootLVItr; 201 G4LogicalVolume* aLV = *rootLVItr;
195 G4Region* aR = *rItr; 202 G4Region* aR = *rItr;
196 203
197 ScanAndSetCouple(aLV,aCouple,aR); 204 ScanAndSetCouple(aLV,aCouple,aR);
198 205
199 // Proceed to the next root logical 206 // Proceed to the next root logical volume in this region
200 ++rootLVItr; 207 ++rootLVItr;
201 } 208 }
202 209
203 // Proceed to next material in this re 210 // Proceed to next material in this region
204 ++mItr; 211 ++mItr;
205 } 212 }
206 } 213 }
207 } 214 }
208 215
209 // Check if sizes of Range/Energy cuts table 216 // Check if sizes of Range/Energy cuts tables are equal to the size of
210 // the couple table. If new couples are made 217 // the couple table. If new couples are made during the previous procedure,
211 // nCouple becomes larger then nTable 218 // nCouple becomes larger then nTable
212 219
213 std::size_t nCouple = coupleTable.size(); 220 std::size_t nCouple = coupleTable.size();
214 std::size_t nTable = energyCutTable[0]->size 221 std::size_t nTable = energyCutTable[0]->size();
215 G4bool newCoupleAppears = nCouple>nTable; 222 G4bool newCoupleAppears = nCouple>nTable;
216 if(newCoupleAppears) 223 if(newCoupleAppears)
217 { 224 {
218 for(std::size_t n=nCouple-nTable; n>0; --n 225 for(std::size_t n=nCouple-nTable; n>0; --n)
219 { 226 {
220 for(std::size_t nn=0; nn< NumberOfG4CutI 227 for(std::size_t nn=0; nn< NumberOfG4CutIndex; ++nn)
221 { 228 {
222 rangeCutTable[nn]->push_back(-1.); 229 rangeCutTable[nn]->push_back(-1.);
223 energyCutTable[nn]->push_back(-1.); 230 energyCutTable[nn]->push_back(-1.);
224 } 231 }
225 } 232 }
226 } 233 }
227 234
228 // resize Range/Energy cuts double vectors i <<
229 if(newCoupleAppears) <<
230 { <<
231 for(std::size_t ix=0; ix<NumberOfG4CutInde <<
232 { <<
233 G4double* rangeVOld = rangeDoubleVector[ <<
234 G4double* energyVOld = energyDoubleVecto <<
235 if(rangeVOld) delete [] rangeVOld; <<
236 if(energyVOld) delete [] energyVOld; <<
237 rangeDoubleVector[ix] = new G4double[(*( <<
238 energyDoubleVector[ix] = new G4double[(* <<
239 } <<
240 } <<
241 } <<
242 <<
243 // ------------------------------------------- <<
244 void G4ProductionCutsTable::UpdateCoupleTable( <<
245 { <<
246 CreateCoupleTables(); <<
247 <<
248 if(firstUse) <<
249 { <<
250 if(G4ParticleTable::GetParticleTable()->Fi <<
251 { <<
252 converters[0] = new G4RToEConvForGamma() <<
253 converters[0]->SetVerboseLevel(GetVerbos <<
254 } <<
255 if(G4ParticleTable::GetParticleTable()->Fi <<
256 { <<
257 converters[1] = new G4RToEConvForElectro <<
258 converters[1]->SetVerboseLevel(GetVerbos <<
259 } <<
260 if(G4ParticleTable::GetParticleTable()->Fi <<
261 { <<
262 converters[2] = new G4RToEConvForPositro <<
263 converters[2]->SetVerboseLevel(GetVerbos <<
264 } <<
265 if(G4ParticleTable::GetParticleTable()->Fi <<
266 { <<
267 converters[3] = new G4RToEConvForProton( <<
268 converters[3]->SetVerboseLevel(GetVerbos <<
269 } <<
270 firstUse = false; <<
271 } <<
272 <<
273 // Update RangeEnergy cuts tables 235 // Update RangeEnergy cuts tables
274 std::size_t idx = 0; 236 std::size_t idx = 0;
275 G4Timer timer; 237 G4Timer timer;
276 if (verboseLevel>2) 238 if (verboseLevel>2)
277 { 239 {
278 timer.Start(); 240 timer.Start();
279 } 241 }
280 for(auto cItr=coupleTable.cbegin(); cItr!=co 242 for(auto cItr=coupleTable.cbegin(); cItr!=coupleTable.cend(); ++cItr)
281 { 243 {
282 G4ProductionCuts* aCut = (*cItr)->GetProdu 244 G4ProductionCuts* aCut = (*cItr)->GetProductionCuts();
283 const G4Material* aMat = (*cItr)->GetMater 245 const G4Material* aMat = (*cItr)->GetMaterial();
284 if((*cItr)->IsRecalcNeeded()) 246 if((*cItr)->IsRecalcNeeded())
285 { 247 {
286 for(std::size_t ptcl=0; ptcl< NumberOfG4 248 for(std::size_t ptcl=0; ptcl< NumberOfG4CutIndex; ++ptcl)
287 { 249 {
288 G4double rCut = aCut->GetProductionCut << 250 G4double rCut = aCut->GetProductionCut(ptcl);
289 (*(rangeCutTable[ptcl]))[idx] = rCut; 251 (*(rangeCutTable[ptcl]))[idx] = rCut;
290 << 252 // if(converters[ptcl] && (*cItr)->IsUsed())
291 if(nullptr != converters[ptcl]) << 253 if(converters[ptcl])
292 { 254 {
293 // check user defined cut in energy << 255 (*(energyCutTable[ptcl]))[idx] = converters[ptcl]->Convert(rCut,aMat);
294 if(nullptr == userEnergyCuts[ptcl] | <<
295 { <<
296 (*(energyCutTable[ptcl]))[idx] = c <<
297 } <<
298 else <<
299 { <<
300 (*(energyCutTable[ptcl]))[idx] = ( <<
301 } <<
302 } 256 }
303 else 257 else
304 { 258 {
305 (*(energyCutTable[ptcl]))[idx] = -1. 259 (*(energyCutTable[ptcl]))[idx] = -1.;
306 } 260 }
307 } 261 }
308 } 262 }
309 ++idx; 263 ++idx;
310 } 264 }
311 if (verboseLevel>2) 265 if (verboseLevel>2)
312 { 266 {
313 timer.Stop(); 267 timer.Stop();
314 G4cout << "G4ProductionCutsTable::UpdateCo 268 G4cout << "G4ProductionCutsTable::UpdateCoupleTable() - "
315 << "Elapsed time for calculation of 269 << "Elapsed time for calculation of energy cuts: " << G4endl;
316 G4cout << timer << G4endl; 270 G4cout << timer << G4endl;
317 } 271 }
318 272
>> 273 // resize Range/Energy cuts double vectors if new couple is made
>> 274 if(newCoupleAppears)
>> 275 {
>> 276 for(std::size_t ix=0; ix<NumberOfG4CutIndex; ++ix)
>> 277 {
>> 278 G4double* rangeVOld = rangeDoubleVector[ix];
>> 279 G4double* energyVOld = energyDoubleVector[ix];
>> 280 if(rangeVOld) delete [] rangeVOld;
>> 281 if(energyVOld) delete [] energyVOld;
>> 282 rangeDoubleVector[ix] = new G4double[(*(rangeCutTable[ix])).size()];
>> 283 energyDoubleVector[ix] = new G4double[(*(energyCutTable[ix])).size()];
>> 284 }
>> 285 }
>> 286
319 // Update Range/Energy cuts double vectors 287 // Update Range/Energy cuts double vectors
320 for(std::size_t ix=0; ix<NumberOfG4CutIndex; 288 for(std::size_t ix=0; ix<NumberOfG4CutIndex; ++ix)
321 { 289 {
322 for(std::size_t ixx=0; ixx<(*(rangeCutTabl 290 for(std::size_t ixx=0; ixx<(*(rangeCutTable[ix])).size(); ++ixx)
323 { 291 {
324 rangeDoubleVector[ix][ixx] = (*(rangeCut 292 rangeDoubleVector[ix][ixx] = (*(rangeCutTable[ix]))[ixx];
325 energyDoubleVector[ix][ixx] = (*(energyC 293 energyDoubleVector[ix][ixx] = (*(energyCutTable[ix]))[ixx];
326 } 294 }
327 } 295 }
328 } 296 }
329 297
330 // ------------------------------------------- 298 // --------------------------------------------------------------------
331 G4double G4ProductionCutsTable::ConvertRangeTo 299 G4double G4ProductionCutsTable::ConvertRangeToEnergy(
332 const G4Part 300 const G4ParticleDefinition* particle,
333 const G4Mate 301 const G4Material* material,
334 G4double 302 G4double range )
335 { 303 {
336 // This method gives energy corresponding to 304 // This method gives energy corresponding to range value
337 305
338 // protection against premature call 306 // protection against premature call
339 if(firstUse) 307 if(firstUse)
340 { 308 {
341 #ifdef G4VERBOSE 309 #ifdef G4VERBOSE
342 if(verboseLevel>0) 310 if(verboseLevel>0)
343 { 311 {
344 G4ExceptionDescription ed; 312 G4ExceptionDescription ed;
345 ed << "Invoked prematurely before it is 313 ed << "Invoked prematurely before it is fully initialized.";
346 G4Exception("G4ProductionCutsTable::Conv 314 G4Exception("G4ProductionCutsTable::ConvertRangeToEnergy()",
347 "CUTS0100", JustWarning, ed) 315 "CUTS0100", JustWarning, ed);
348 } 316 }
349 #endif 317 #endif
350 return -1.0; 318 return -1.0;
351 } 319 }
352 320
353 // check material 321 // check material
354 if (material == nullptr) return -1.0; 322 if (material == nullptr) return -1.0;
355 323
356 // check range 324 // check range
357 if (range == 0.0) return 0.0; 325 if (range == 0.0) return 0.0;
358 if (range <0.0) return -1.0; 326 if (range <0.0) return -1.0;
359 327
360 // check particle 328 // check particle
361 G4int index = G4ProductionCuts::GetIndex(par 329 G4int index = G4ProductionCuts::GetIndex(particle);
362 330
363 if (index<0 || converters[index] == nullptr) 331 if (index<0 || converters[index] == nullptr)
364 { 332 {
365 #ifdef G4VERBOSE 333 #ifdef G4VERBOSE
366 if(verboseLevel>0) 334 if(verboseLevel>0)
367 { 335 {
368 G4ExceptionDescription ed; 336 G4ExceptionDescription ed;
369 ed << "Invoked "; 337 ed << "Invoked ";
370 if(particle != nullptr) 338 if(particle != nullptr)
371 { ed << "for particle <" << particle->Ge 339 { ed << "for particle <" << particle->GetParticleName() << ">."; }
372 else 340 else
373 { ed << "without valid particle pointer. 341 { ed << "without valid particle pointer."; }
374 G4Exception("G4ProductionCutsTable::Conv 342 G4Exception("G4ProductionCutsTable::ConvertRangeToEnergy()",
375 "CUTS0101", JustWarning, ed) 343 "CUTS0101", JustWarning, ed);
376 } 344 }
377 #endif 345 #endif
378 return -1.0; 346 return -1.0;
379 } 347 }
380 348
381 return converters[index]->Convert(range, mat 349 return converters[index]->Convert(range, material);
382 } 350 }
383 351
384 // ------------------------------------------- 352 // --------------------------------------------------------------------
385 void G4ProductionCutsTable::ResetConverters() 353 void G4ProductionCutsTable::ResetConverters()
386 {} 354 {}
387 355
388 // ------------------------------------------- 356 // --------------------------------------------------------------------
389 void G4ProductionCutsTable::SetEnergyRange(G4d 357 void G4ProductionCutsTable::SetEnergyRange(G4double lowedge, G4double highedge)
390 { 358 {
391 G4VRangeToEnergyConverter::SetEnergyRange(lo 359 G4VRangeToEnergyConverter::SetEnergyRange(lowedge,highedge);
392 } 360 }
393 361
394 // ------------------------------------------- 362 // --------------------------------------------------------------------
395 G4double G4ProductionCutsTable::GetLowEdgeEne 363 G4double G4ProductionCutsTable::GetLowEdgeEnergy() const
396 { 364 {
397 return G4VRangeToEnergyConverter::GetLowEdge 365 return G4VRangeToEnergyConverter::GetLowEdgeEnergy();
398 } 366 }
399 367
400 // ------------------------------------------- 368 // --------------------------------------------------------------------
401 G4double G4ProductionCutsTable::GetHighEdgeEne 369 G4double G4ProductionCutsTable::GetHighEdgeEnergy() const
402 { 370 {
403 return G4VRangeToEnergyConverter::GetHighEdg 371 return G4VRangeToEnergyConverter::GetHighEdgeEnergy();
404 } 372 }
405 373
406 // ------------------------------------------- 374 // --------------------------------------------------------------------
407 void G4ProductionCutsTable::ScanAndSetCouple(G 375 void G4ProductionCutsTable::ScanAndSetCouple(G4LogicalVolume* aLV,
408 G 376 G4MaterialCutsCouple* aCouple,
409 G 377 G4Region* aRegion)
410 { 378 {
411 // Check whether or not this logical volume 379 // Check whether or not this logical volume belongs to the same region
412 if((aRegion!=nullptr) && aLV->GetRegion()!=a 380 if((aRegion!=nullptr) && aLV->GetRegion()!=aRegion) return;
413 381
414 // Check if this particular volume has a mat 382 // Check if this particular volume has a material matched to the couple
415 if(aLV->GetMaterial()==aCouple->GetMaterial( 383 if(aLV->GetMaterial()==aCouple->GetMaterial())
416 { 384 {
417 aLV->SetMaterialCutsCouple(aCouple); 385 aLV->SetMaterialCutsCouple(aCouple);
418 } 386 }
419 387
420 std::size_t noDaughters = aLV->GetNoDaughter 388 std::size_t noDaughters = aLV->GetNoDaughters();
421 if(noDaughters==0) return; 389 if(noDaughters==0) return;
422 390
423 // Loop over daughters with same region 391 // Loop over daughters with same region
424 for(std::size_t i=0; i<noDaughters; ++i) 392 for(std::size_t i=0; i<noDaughters; ++i)
425 { 393 {
426 G4LogicalVolume* daughterLVol = aLV->GetDa 394 G4LogicalVolume* daughterLVol = aLV->GetDaughter(i)->GetLogicalVolume();
427 ScanAndSetCouple(daughterLVol,aCouple,aReg 395 ScanAndSetCouple(daughterLVol,aCouple,aRegion);
428 } 396 }
429 } 397 }
430 398
431 // ------------------------------------------- 399 // --------------------------------------------------------------------
432 void G4ProductionCutsTable::DumpCouples() cons 400 void G4ProductionCutsTable::DumpCouples() const
433 { 401 {
434 G4cout << G4endl; 402 G4cout << G4endl;
435 G4cout << "========= Table of registered cou 403 G4cout << "========= Table of registered couples ============================"
436 << G4endl; 404 << G4endl;
437 for(auto cItr=coupleTable.cbegin(); cItr!=co 405 for(auto cItr=coupleTable.cbegin(); cItr!=coupleTable.cend(); ++cItr)
438 { 406 {
439 G4MaterialCutsCouple* aCouple = (*cItr); 407 G4MaterialCutsCouple* aCouple = (*cItr);
440 G4ProductionCuts* aCut = aCouple->GetProdu 408 G4ProductionCuts* aCut = aCouple->GetProductionCuts();
441 G4cout << G4endl; 409 G4cout << G4endl;
442 G4cout << "Index : " << aCouple->GetIndex( 410 G4cout << "Index : " << aCouple->GetIndex()
443 << " used in the geometry : "; 411 << " used in the geometry : ";
444 if(aCouple->IsUsed()) G4cout << "Yes"; 412 if(aCouple->IsUsed()) G4cout << "Yes";
445 else G4cout << "No "; 413 else G4cout << "No ";
446 //// G4cout << " recalculation needed : 414 //// G4cout << " recalculation needed : ";
447 //// if(aCouple->IsRecalcNeeded()) G4cout < 415 //// if(aCouple->IsRecalcNeeded()) G4cout << "Yes";
448 //// else G4cout < 416 //// else G4cout << "No ";
449 G4cout << G4endl; 417 G4cout << G4endl;
450 G4cout << " Material : " << aCouple->GetMa 418 G4cout << " Material : " << aCouple->GetMaterial()->GetName() << G4endl;
451 G4cout << " Range cuts : " 419 G4cout << " Range cuts : "
452 << " gamma " << G4BestUnit(aCut->G 420 << " gamma " << G4BestUnit(aCut->GetProductionCut("gamma"),"Length")
453 << " e- " << G4BestUnit(aCut->G 421 << " e- " << G4BestUnit(aCut->GetProductionCut("e-"),"Length")
454 << " e+ " << G4BestUnit(aCut->G 422 << " e+ " << G4BestUnit(aCut->GetProductionCut("e+"),"Length")
455 << " proton " << G4BestUnit(aCut->G 423 << " proton " << G4BestUnit(aCut->GetProductionCut("proton"),"Length")
456 << G4endl; 424 << G4endl;
457 G4cout << " Energy thresholds : " ; 425 G4cout << " Energy thresholds : " ;
458 //// if(aCouple->IsRecalcNeeded()) { 426 //// if(aCouple->IsRecalcNeeded()) {
459 //// G4cout << " is not ready to print"; 427 //// G4cout << " is not ready to print";
460 //// } else { 428 //// } else {
461 G4cout << " gamma " << G4BestUnit((*(ener 429 G4cout << " gamma " << G4BestUnit((*(energyCutTable[0]))[aCouple->GetIndex()],"Energy")
462 << " e- " << G4BestUnit((*(ener 430 << " e- " << G4BestUnit((*(energyCutTable[1]))[aCouple->GetIndex()],"Energy")
463 << " e+ " << G4BestUnit((*(ener 431 << " e+ " << G4BestUnit((*(energyCutTable[2]))[aCouple->GetIndex()],"Energy")
464 << " proton " << G4BestUnit((*(ener 432 << " proton " << G4BestUnit((*(energyCutTable[3]))[aCouple->GetIndex()],"Energy");
465 //// } 433 //// }
466 G4cout << G4endl; 434 G4cout << G4endl;
467 435
468 if(aCouple->IsUsed()) 436 if(aCouple->IsUsed())
469 { 437 {
470 G4cout << " Region(s) which use this cou 438 G4cout << " Region(s) which use this couple : " << G4endl;
471 for(auto rItr=fG4RegionStore->cbegin(); 439 for(auto rItr=fG4RegionStore->cbegin();
472 rItr!=fG4RegionStore->cend(); + 440 rItr!=fG4RegionStore->cend(); ++rItr)
473 { 441 {
474 if (IsCoupleUsedInTheRegion(aCouple, * 442 if (IsCoupleUsedInTheRegion(aCouple, *rItr) )
475 { 443 {
476 G4cout << " " << (*rItr)->GetName 444 G4cout << " " << (*rItr)->GetName() << G4endl;
477 } 445 }
478 } 446 }
479 } 447 }
480 } 448 }
481 G4cout << G4endl; 449 G4cout << G4endl;
482 G4cout << "================================= 450 G4cout << "==================================================================" << G4endl;
483 G4cout << G4endl; 451 G4cout << G4endl;
484 } 452 }
485 453
486 // ------------------------------------------- 454 // --------------------------------------------------------------------
487 G4bool G4ProductionCutsTable::StoreCutsTable( 455 G4bool G4ProductionCutsTable::StoreCutsTable(const G4String& dir,
488 456 G4bool ascii)
489 { 457 {
490 // Store cuts and material information in fi 458 // Store cuts and material information in files under the specified directory
491 459
492 if (!StoreMaterialInfo(dir, ascii)) return f 460 if (!StoreMaterialInfo(dir, ascii)) return false;
493 if (!StoreMaterialCutsCoupleInfo(dir, ascii) 461 if (!StoreMaterialCutsCoupleInfo(dir, ascii)) return false;
494 if (!StoreCutsInfo(dir, ascii)) return false 462 if (!StoreCutsInfo(dir, ascii)) return false;
495 463
496 #ifdef G4VERBOSE 464 #ifdef G4VERBOSE
497 if (verboseLevel >2) 465 if (verboseLevel >2)
498 { 466 {
499 G4cout << "G4ProductionCutsTable::StoreCut 467 G4cout << "G4ProductionCutsTable::StoreCutsTable()" << G4endl;
500 G4cout << " Material/Cuts information have 468 G4cout << " Material/Cuts information have been successfully stored ";
501 if (ascii) 469 if (ascii)
502 { 470 {
503 G4cout << " in Ascii mode "; 471 G4cout << " in Ascii mode ";
504 } 472 }
505 else 473 else
506 { 474 {
507 G4cout << " in Binary mode "; 475 G4cout << " in Binary mode ";
508 } 476 }
509 G4cout << " under " << dir << G4endl; 477 G4cout << " under " << dir << G4endl;
510 } 478 }
511 #endif 479 #endif
512 return true; 480 return true;
513 } 481 }
514 482
515 // ------------------------------------------- 483 // --------------------------------------------------------------------
516 G4bool G4ProductionCutsTable::RetrieveCutsTab 484 G4bool G4ProductionCutsTable::RetrieveCutsTable(const G4String& dir,
517 485 G4bool ascii)
518 { 486 {
519 if (!CheckForRetrieveCutsTable(dir, ascii)) 487 if (!CheckForRetrieveCutsTable(dir, ascii)) return false;
520 if (!RetrieveCutsInfo(dir, ascii)) return fa 488 if (!RetrieveCutsInfo(dir, ascii)) return false;
521 #ifdef G4VERBOSE 489 #ifdef G4VERBOSE
522 if (verboseLevel >2) 490 if (verboseLevel >2)
523 { 491 {
524 G4cout << "G4ProductionCutsTable::Retrieve 492 G4cout << "G4ProductionCutsTable::RetrieveCutsTable()" << G4endl;
525 G4cout << " Material/Cuts information have 493 G4cout << " Material/Cuts information have been successfully retrieved ";
526 if (ascii) 494 if (ascii)
527 { 495 {
528 G4cout << " in Ascii mode "; 496 G4cout << " in Ascii mode ";
529 } 497 }
530 else 498 else
531 { 499 {
532 G4cout << " in Binary mode "; 500 G4cout << " in Binary mode ";
533 } 501 }
534 G4cout << " under " << dir << G4endl; 502 G4cout << " under " << dir << G4endl;
535 } 503 }
536 #endif 504 #endif
537 return true; 505 return true;
538 } 506 }
539 507
540 // ------------------------------------------- 508 // --------------------------------------------------------------------
541 G4bool 509 G4bool
542 G4ProductionCutsTable::CheckForRetrieveCutsTab 510 G4ProductionCutsTable::CheckForRetrieveCutsTable(const G4String& directory,
543 511 G4bool ascii)
544 { 512 {
545 // check stored material and cut values are 513 // check stored material and cut values are consistent
546 // with the current detector setup 514 // with the current detector setup
547 515
548 G4cerr << "G4ProductionCutsTable::CheckForRe 516 G4cerr << "G4ProductionCutsTable::CheckForRetrieveCutsTable()"<< G4endl;
549 // isNeedForRestoreCoupleInfo = false; 517 // isNeedForRestoreCoupleInfo = false;
550 if (!CheckMaterialInfo(directory, ascii)) re 518 if (!CheckMaterialInfo(directory, ascii)) return false;
551 if (verboseLevel >2) 519 if (verboseLevel >2)
552 { 520 {
553 G4cerr << "G4ProductionCutsTable::CheckM 521 G4cerr << "G4ProductionCutsTable::CheckMaterialInfo passed !!"<< G4endl;
554 } 522 }
555 if (!CheckMaterialCutsCoupleInfo(directory, 523 if (!CheckMaterialCutsCoupleInfo(directory, ascii)) return false;
556 if (verboseLevel >2) 524 if (verboseLevel >2)
557 { 525 {
558 G4cerr << "G4ProductionCutsTable::CheckMat 526 G4cerr << "G4ProductionCutsTable::CheckMaterialCutsCoupleInfo passed !!"
559 << G4endl; 527 << G4endl;
560 } 528 }
561 return true; 529 return true;
562 } 530 }
563 531
564 // ------------------------------------------- 532 // --------------------------------------------------------------------
565 G4bool G4ProductionCutsTable::StoreMaterialInf 533 G4bool G4ProductionCutsTable::StoreMaterialInfo(const G4String& directory,
566 534 G4bool ascii)
567 { 535 {
568 // Store material information in files under 536 // Store material information in files under the specified directory
569 537
570 const G4String fileName = directory + "/" + 538 const G4String fileName = directory + "/" + "material.dat";
571 const G4String key = "MATERIAL-V3.0"; 539 const G4String key = "MATERIAL-V3.0";
572 std::ofstream fOut; 540 std::ofstream fOut;
573 541
574 // open output file 542 // open output file
575 if (!ascii ) fOut.open(fileName,std::ios::o 543 if (!ascii ) fOut.open(fileName,std::ios::out|std::ios::binary);
576 else fOut.open(fileName,std::ios::o 544 else fOut.open(fileName,std::ios::out);
577 545
578 // check if the file has been opened success 546 // check if the file has been opened successfully
579 if (!fOut) 547 if (!fOut)
580 { 548 {
581 #ifdef G4VERBOSE 549 #ifdef G4VERBOSE
582 if (verboseLevel>0) 550 if (verboseLevel>0)
583 { 551 {
584 G4cerr << "G4ProductionCutsTable::StoreM 552 G4cerr << "G4ProductionCutsTable::StoreMaterialInfo() - ";
585 G4cerr << "Cannot open file: " << fileNa 553 G4cerr << "Cannot open file: " << fileName << G4endl;
586 } 554 }
587 #endif 555 #endif
588 G4Exception( "G4ProductionCutsTable::Store 556 G4Exception( "G4ProductionCutsTable::StoreMaterialInfo()",
589 "ProcCuts102", JustWarning, " 557 "ProcCuts102", JustWarning, "Cannot open file!");
590 return false; 558 return false;
591 } 559 }
592 560
593 const G4MaterialTable* matTable = G4Material 561 const G4MaterialTable* matTable = G4Material::GetMaterialTable();
594 // number of materials in the table 562 // number of materials in the table
595 G4int numberOfMaterial = (G4int)matTable->si << 563 G4int numberOfMaterial = matTable->size();
596 564
597 if (ascii) 565 if (ascii)
598 { 566 {
599 /////////////// ASCII mode ////////////// 567 /////////////// ASCII mode /////////////////
600 // key word 568 // key word
601 fOut << key << G4endl; 569 fOut << key << G4endl;
602 570
603 // number of materials in the table 571 // number of materials in the table
604 fOut << numberOfMaterial << G4endl; 572 fOut << numberOfMaterial << G4endl;
605 573
606 fOut.setf(std::ios::scientific); 574 fOut.setf(std::ios::scientific);
607 575
608 // material name and density 576 // material name and density
609 for (std::size_t idx=0; static_cast<G4int> 577 for (std::size_t idx=0; static_cast<G4int>(idx)<numberOfMaterial; ++idx)
610 { 578 {
611 fOut << std::setw(FixedStringLengthForSt 579 fOut << std::setw(FixedStringLengthForStore)
612 << ((*matTable)[idx])->GetName(); 580 << ((*matTable)[idx])->GetName();
613 fOut << std::setw(FixedStringLengthForSt 581 fOut << std::setw(FixedStringLengthForStore)
614 << ((*matTable)[idx])->GetDensity() 582 << ((*matTable)[idx])->GetDensity()/(g/cm3) << G4endl;
615 } 583 }
616 584
617 fOut.unsetf(std::ios::scientific); 585 fOut.unsetf(std::ios::scientific);
618 586
619 } 587 }
620 else 588 else
621 { 589 {
622 /////////////// Binary mode ///////////// 590 /////////////// Binary mode /////////////////
623 char temp[FixedStringLengthForStore]; 591 char temp[FixedStringLengthForStore];
624 std::size_t i; 592 std::size_t i;
625 593
626 // key word 594 // key word
627 for (i=0; i<FixedStringLengthForStore; ++i 595 for (i=0; i<FixedStringLengthForStore; ++i)
628 { 596 {
629 temp[i] = '\0'; 597 temp[i] = '\0';
630 } 598 }
631 for (i=0; i<key.length() && i<FixedStringL 599 for (i=0; i<key.length() && i<FixedStringLengthForStore-1; ++i)
632 { 600 {
633 temp[i]=key[(G4int)i]; << 601 temp[i]=key[i];
634 } 602 }
635 fOut.write(temp, FixedStringLengthForStore 603 fOut.write(temp, FixedStringLengthForStore);
636 604
637 // number of materials in the table 605 // number of materials in the table
638 fOut.write( (char*)(&numberOfMaterial), si 606 fOut.write( (char*)(&numberOfMaterial), sizeof(G4int));
639 607
640 // material name and density 608 // material name and density
641 for (std::size_t imat=0; static_cast<G4int 609 for (std::size_t imat=0; static_cast<G4int>(imat)<numberOfMaterial; ++imat)
642 { 610 {
643 G4String name = ((*matTable)[imat])->Ge 611 G4String name = ((*matTable)[imat])->GetName();
644 G4double density = ((*matTable)[imat])-> 612 G4double density = ((*matTable)[imat])->GetDensity();
645 for (i=0; i<FixedStringLengthForStore; + 613 for (i=0; i<FixedStringLengthForStore; ++i)
646 temp[i] = '\0'; 614 temp[i] = '\0';
647 for (i=0; i<name.length() && i<FixedStri 615 for (i=0; i<name.length() && i<FixedStringLengthForStore-1; ++i)
648 temp[i]=name[(G4int)i]; << 616 temp[i]=name[i];
649 fOut.write(temp, FixedStringLengthForSto 617 fOut.write(temp, FixedStringLengthForStore);
650 fOut.write( (char*)(&density), sizeof(G4 618 fOut.write( (char*)(&density), sizeof(G4double));
651 } 619 }
652 } 620 }
653 621
654 fOut.close(); 622 fOut.close();
655 return true; 623 return true;
656 } 624 }
657 625
658 // ------------------------------------------- 626 // --------------------------------------------------------------------
659 G4bool G4ProductionCutsTable::CheckMaterialInf 627 G4bool G4ProductionCutsTable::CheckMaterialInfo(const G4String& directory,
660 628 G4bool ascii)
661 { 629 {
662 // Check stored material is consistent with 630 // Check stored material is consistent with the current detector setup
663 631
664 const G4String fileName = directory + "/" + 632 const G4String fileName = directory + "/" + "material.dat";
665 const G4String key = "MATERIAL-V3.0"; 633 const G4String key = "MATERIAL-V3.0";
666 std::ifstream fIn; 634 std::ifstream fIn;
667 635
668 // open input file 636 // open input file
669 if (!ascii ) fIn.open(fileName,std::ios::in| 637 if (!ascii ) fIn.open(fileName,std::ios::in|std::ios::binary);
670 else fIn.open(fileName,std::ios::in) 638 else fIn.open(fileName,std::ios::in);
671 639
672 // check if the file has been opened success 640 // check if the file has been opened successfully
673 if (!fIn) 641 if (!fIn)
674 { 642 {
675 #ifdef G4VERBOSE 643 #ifdef G4VERBOSE
676 if (verboseLevel >0) 644 if (verboseLevel >0)
677 { 645 {
678 G4cerr << "G4ProductionCutsTable::CheckM 646 G4cerr << "G4ProductionCutsTable::CheckMaterialInfo() - ";
679 G4cerr << "Cannot open file: " << fileNa 647 G4cerr << "Cannot open file: " << fileName << G4endl;
680 } 648 }
681 #endif 649 #endif
682 G4Exception( "G4ProductionCutsTable::Check 650 G4Exception( "G4ProductionCutsTable::CheckMaterialInfo()",
683 "ProcCuts102", JustWarning, " 651 "ProcCuts102", JustWarning, "Cannot open file!");
684 return false; 652 return false;
685 } 653 }
686 654
687 char temp[FixedStringLengthForStore]; 655 char temp[FixedStringLengthForStore];
688 656
689 // key word 657 // key word
690 G4String keyword; 658 G4String keyword;
691 if (ascii) 659 if (ascii)
692 { 660 {
693 fIn >> keyword; 661 fIn >> keyword;
694 } 662 }
695 else 663 else
696 { 664 {
697 fIn.read(temp, FixedStringLengthForStore); 665 fIn.read(temp, FixedStringLengthForStore);
698 keyword = (const char*)(temp); 666 keyword = (const char*)(temp);
699 } 667 }
700 if (key!=keyword) 668 if (key!=keyword)
701 { 669 {
702 #ifdef G4VERBOSE 670 #ifdef G4VERBOSE
703 if (verboseLevel >0) 671 if (verboseLevel >0)
704 { 672 {
705 G4cerr << "G4ProductionCutsTable::CheckM 673 G4cerr << "G4ProductionCutsTable::CheckMaterialInfo() - ";
706 G4cerr << "Key word in " << fileName << 674 G4cerr << "Key word in " << fileName << "= " << keyword ;
707 G4cerr <<"( should be "<< key << ")" < 675 G4cerr <<"( should be "<< key << ")" <<G4endl;
708 } 676 }
709 #endif 677 #endif
710 G4Exception( "G4ProductionCutsTable::Check 678 G4Exception( "G4ProductionCutsTable::CheckMaterialInfo()",
711 "ProcCuts103", JustWarning, " 679 "ProcCuts103", JustWarning, "Bad Data Format");
712 return false; 680 return false;
713 } 681 }
714 682
715 // number of materials in the table 683 // number of materials in the table
716 G4int nmat; 684 G4int nmat;
717 if (ascii) 685 if (ascii)
718 { 686 {
719 fIn >> nmat; 687 fIn >> nmat;
720 } 688 }
721 else 689 else
722 { 690 {
723 fIn.read( (char*)(&nmat), sizeof(G4int)); 691 fIn.read( (char*)(&nmat), sizeof(G4int));
724 } 692 }
725 if ((nmat<=0) || (nmat >100000)) 693 if ((nmat<=0) || (nmat >100000))
726 { 694 {
727 G4Exception( "G4ProductionCutsTable::Check 695 G4Exception( "G4ProductionCutsTable::CheckMaterialInfo()",
728 "ProcCuts108", JustWarning, 696 "ProcCuts108", JustWarning,
729 "Number of materials is less 697 "Number of materials is less than zero or too big");
730 return false; 698 return false;
731 } 699 }
732 700
733 // list of material 701 // list of material
734 for (G4int idx=0; idx<nmat ; ++idx) 702 for (G4int idx=0; idx<nmat ; ++idx)
735 { 703 {
736 // check eof 704 // check eof
737 if(fIn.eof()) 705 if(fIn.eof())
738 { 706 {
739 #ifdef G4VERBOSE 707 #ifdef G4VERBOSE
740 if (verboseLevel >0) 708 if (verboseLevel >0)
741 { 709 {
742 G4cout << "G4ProductionCutsTable::Chec 710 G4cout << "G4ProductionCutsTable::CheckMaterialInfo() - ";
743 G4cout << "Encountered End of File " ; 711 G4cout << "Encountered End of File " ;
744 G4cout << " at " << idx+1 << "th mate 712 G4cout << " at " << idx+1 << "th material "<< G4endl;
745 } 713 }
746 #endif 714 #endif
747 fIn.close(); 715 fIn.close();
748 return false; 716 return false;
749 } 717 }
750 718
751 // check material name and density 719 // check material name and density
752 char name[FixedStringLengthForStore]; 720 char name[FixedStringLengthForStore];
753 G4double density; 721 G4double density;
754 if (ascii) 722 if (ascii)
755 { 723 {
756 fIn >> name >> density; 724 fIn >> name >> density;
757 density *= (g/cm3); 725 density *= (g/cm3);
758 726
759 } 727 }
760 else 728 else
761 { 729 {
762 fIn.read(name, FixedStringLengthForStore 730 fIn.read(name, FixedStringLengthForStore);
763 fIn.read((char*)(&density), sizeof(G4dou 731 fIn.read((char*)(&density), sizeof(G4double));
764 } 732 }
765 if (fIn.fail()) 733 if (fIn.fail())
766 { 734 {
767 #ifdef G4VERBOSE 735 #ifdef G4VERBOSE
768 if (verboseLevel >0) 736 if (verboseLevel >0)
769 { 737 {
770 G4cerr << "G4ProductionCutsTable::Chec 738 G4cerr << "G4ProductionCutsTable::CheckMaterialInfo() - ";
771 G4cerr << "Bad data format "; 739 G4cerr << "Bad data format ";
772 G4cerr << " at " << idx+1 << "th mate 740 G4cerr << " at " << idx+1 << "th material "<< G4endl;
773 } 741 }
774 #endif 742 #endif
775 G4Exception( "G4ProductionCutsTable::Che 743 G4Exception( "G4ProductionCutsTable::CheckMaterialInfo()",
776 "ProcCuts103", JustWarning, 744 "ProcCuts103", JustWarning, "Bad Data Format");
777 fIn.close(); 745 fIn.close();
778 return false; 746 return false;
779 } 747 }
780 748
781 G4Material* aMaterial = G4Material::GetMat 749 G4Material* aMaterial = G4Material::GetMaterial(name);
782 if (aMaterial == nullptr ) continue; 750 if (aMaterial == nullptr ) continue;
783 751
784 G4double ratio = std::fabs(density/aMateri 752 G4double ratio = std::fabs(density/aMaterial->GetDensity() );
785 if ((0.999>ratio) || (ratio>1.001) ) 753 if ((0.999>ratio) || (ratio>1.001) )
786 { 754 {
787 #ifdef G4VERBOSE 755 #ifdef G4VERBOSE
788 if (verboseLevel >0) 756 if (verboseLevel >0)
789 { 757 {
790 G4cerr << "G4ProductionCutsTable::Chec 758 G4cerr << "G4ProductionCutsTable::CheckMaterialInfo() - ";
791 G4cerr << "Inconsistent material densi 759 G4cerr << "Inconsistent material density" << G4endl;;
792 G4cerr << " at " << idx+1 << "th mate 760 G4cerr << " at " << idx+1 << "th material "<< G4endl;
793 G4cerr << "Name: " << name << G4endl 761 G4cerr << "Name: " << name << G4endl;
794 G4cerr << "Density:" << std::setiosfla 762 G4cerr << "Density:" << std::setiosflags(std::ios::scientific)
795 << density / (g/cm3) ; 763 << density / (g/cm3) ;
796 G4cerr << "(should be " << aMaterial-> 764 G4cerr << "(should be " << aMaterial->GetDensity()/(g/cm3)<< ")"
797 << " [g/cm3]"<< G4endl; 765 << " [g/cm3]"<< G4endl;
798 G4cerr << std::resetiosflags(std::ios: 766 G4cerr << std::resetiosflags(std::ios::scientific);
799 } 767 }
800 #endif 768 #endif
801 G4Exception( "G4ProductionCutsTable::Che 769 G4Exception( "G4ProductionCutsTable::CheckMaterialInfo()",
802 "ProcCuts104", JustWarning, 770 "ProcCuts104", JustWarning, "Inconsistent material density");
803 fIn.close(); 771 fIn.close();
804 return false; 772 return false;
805 } 773 }
806 } 774 }
807 775
808 fIn.close(); 776 fIn.close();
809 return true; 777 return true;
810 } 778 }
811 779
812 // ------------------------------------------- 780 // --------------------------------------------------------------------
813 G4bool 781 G4bool
814 G4ProductionCutsTable::StoreMaterialCutsCouple 782 G4ProductionCutsTable::StoreMaterialCutsCoupleInfo(const G4String& directory,
815 783 G4bool ascii)
816 { 784 {
817 // Store materialCutsCouple information in f 785 // Store materialCutsCouple information in files under the specified directory
818 786
819 const G4String fileName = directory + "/" + 787 const G4String fileName = directory + "/" + "couple.dat";
820 const G4String key = "COUPLE-V3.0"; 788 const G4String key = "COUPLE-V3.0";
821 std::ofstream fOut; 789 std::ofstream fOut;
822 char temp[FixedStringLengthForStore]; 790 char temp[FixedStringLengthForStore];
823 791
824 // open output file 792 // open output file
825 if (!ascii ) fOut.open(fileName,std::ios::ou 793 if (!ascii ) fOut.open(fileName,std::ios::out|std::ios::binary);
826 else fOut.open(fileName,std::ios::ou 794 else fOut.open(fileName,std::ios::out);
827 795
828 796
829 // check if the file has been opened success 797 // check if the file has been opened successfully
830 if (!fOut) 798 if (!fOut)
831 { 799 {
832 #ifdef G4VERBOSE 800 #ifdef G4VERBOSE
833 if (verboseLevel >0) 801 if (verboseLevel >0)
834 { 802 {
835 G4cerr << "G4ProductionCutsTable::StoreM 803 G4cerr << "G4ProductionCutsTable::StoreMaterialCutsCoupleInfo() - ";
836 G4cerr << "Cannot open file: " << fileNa 804 G4cerr << "Cannot open file: " << fileName << G4endl;
837 } 805 }
838 #endif 806 #endif
839 G4Exception( "G4ProductionCutsTable::Store 807 G4Exception( "G4ProductionCutsTable::StoreMaterialCutsCoupleInfo()",
840 "ProcCuts102", 808 "ProcCuts102",
841 JustWarning, "Cannot open fil 809 JustWarning, "Cannot open file!");
842 return false; 810 return false;
843 } 811 }
844 G4int numberOfCouples = (G4int)coupleTable.s << 812 G4int numberOfCouples = coupleTable.size();
845 if (ascii) 813 if (ascii)
846 { 814 {
847 /////////////// ASCII mode ////////////// 815 /////////////// ASCII mode /////////////////
848 // key word 816 // key word
849 fOut << std::setw(FixedStringLengthForStor 817 fOut << std::setw(FixedStringLengthForStore) << key << G4endl;
850 818
851 // number of couples in the table 819 // number of couples in the table
852 fOut << numberOfCouples << G4endl; 820 fOut << numberOfCouples << G4endl;
853 } 821 }
854 else 822 else
855 { 823 {
856 /////////////// Binary mode ///////////// 824 /////////////// Binary mode /////////////////
857 // key word 825 // key word
858 std::size_t i; 826 std::size_t i;
859 for (i=0; i<FixedStringLengthForStore; ++i 827 for (i=0; i<FixedStringLengthForStore; ++i)
860 temp[i] = '\0'; 828 temp[i] = '\0';
861 for (i=0; i<key.length() && i<FixedStringL 829 for (i=0; i<key.length() && i<FixedStringLengthForStore-1; ++i)
862 temp[i]=key[(G4int)i]; << 830 temp[i]=key[i];
863 fOut.write(temp, FixedStringLengthForStore 831 fOut.write(temp, FixedStringLengthForStore);
864 832
865 // number of couples in the table 833 // number of couples in the table
866 fOut.write( (char*)(&numberOfCouples), siz 834 fOut.write( (char*)(&numberOfCouples), sizeof(G4int));
867 } 835 }
868 836
869 // Loop over all couples 837 // Loop over all couples
870 for (auto cItr=coupleTable.cbegin(); cItr!=c 838 for (auto cItr=coupleTable.cbegin(); cItr!=coupleTable.cend(); ++cItr)
871 { 839 {
872 G4MaterialCutsCouple* aCouple = (*cItr); 840 G4MaterialCutsCouple* aCouple = (*cItr);
873 G4int index = aCouple->GetIndex(); 841 G4int index = aCouple->GetIndex();
874 // cut value 842 // cut value
875 G4ProductionCuts* aCut = aCouple->GetProdu 843 G4ProductionCuts* aCut = aCouple->GetProductionCuts();
876 G4double cutValues[NumberOfG4CutIndex]; 844 G4double cutValues[NumberOfG4CutIndex];
877 for (std::size_t idx=0; idx <NumberOfG4Cut 845 for (std::size_t idx=0; idx <NumberOfG4CutIndex; ++idx)
878 { 846 {
879 cutValues[idx] = aCut->GetProductionCut( << 847 cutValues[idx] = aCut->GetProductionCut(idx);
880 } 848 }
881 // material/region info 849 // material/region info
882 G4String materialName = aCouple->GetMateri 850 G4String materialName = aCouple->GetMaterial()->GetName();
883 G4String regionName = "NONE"; 851 G4String regionName = "NONE";
884 if (aCouple->IsUsed()) 852 if (aCouple->IsUsed())
885 { 853 {
886 for(auto rItr=fG4RegionStore->cbegin(); 854 for(auto rItr=fG4RegionStore->cbegin();
887 rItr!=fG4RegionStore->cend(); + 855 rItr!=fG4RegionStore->cend(); ++rItr)
888 { 856 {
889 if (IsCoupleUsedInTheRegion(aCouple, * 857 if (IsCoupleUsedInTheRegion(aCouple, *rItr))
890 { 858 {
891 regionName = (*rItr)->GetName(); 859 regionName = (*rItr)->GetName();
892 break; 860 break;
893 } 861 }
894 } 862 }
895 } 863 }
896 864
897 if (ascii) 865 if (ascii)
898 { 866 {
899 /////////////// ASCII mode //////////// 867 /////////////// ASCII mode /////////////////
900 // index number 868 // index number
901 fOut << index << G4endl; 869 fOut << index << G4endl;
902 870
903 // material name 871 // material name
904 fOut << std::setw(FixedStringLengthForSt 872 fOut << std::setw(FixedStringLengthForStore) << materialName<< G4endl;
905 873
906 // region name 874 // region name
907 fOut << std::setw(FixedStringLengthForSt 875 fOut << std::setw(FixedStringLengthForStore) << regionName<< G4endl;
908 876
909 fOut.setf(std::ios::scientific); 877 fOut.setf(std::ios::scientific);
910 // cut values 878 // cut values
911 for (std::size_t idx=0; idx< NumberOfG4C 879 for (std::size_t idx=0; idx< NumberOfG4CutIndex; ++idx)
912 { 880 {
913 fOut << std::setw(FixedStringLengthFor 881 fOut << std::setw(FixedStringLengthForStore) << cutValues[idx]/(mm)
914 << G4endl; 882 << G4endl;
915 } 883 }
916 fOut.unsetf(std::ios::scientific); 884 fOut.unsetf(std::ios::scientific);
917 885
918 } 886 }
919 else 887 else
920 { 888 {
921 /////////////// Binary mode /////////// 889 /////////////// Binary mode /////////////////
922 // index 890 // index
923 fOut.write( (char*)(&index), sizeof(G4in 891 fOut.write( (char*)(&index), sizeof(G4int));
924 892
925 // material name 893 // material name
926 std::size_t i; 894 std::size_t i;
927 for (i=0; i<FixedStringLengthForStore; + 895 for (i=0; i<FixedStringLengthForStore; ++i)
928 temp[i] = '\0'; 896 temp[i] = '\0';
929 for (i=0; i<materialName.length() && i<F 897 for (i=0; i<materialName.length() && i<FixedStringLengthForStore-1; ++i)
930 temp[i]=materialName[(G4int)i]; << 898 temp[i]=materialName[i];
931 fOut.write(temp, FixedStringLengthForSto 899 fOut.write(temp, FixedStringLengthForStore);
932 900
933 // region name 901 // region name
934 for (i=0; i<FixedStringLengthForStore; + 902 for (i=0; i<FixedStringLengthForStore; ++i)
935 temp[i] = '\0'; 903 temp[i] = '\0';
936 for (i=0; i<regionName.length() && i<Fix 904 for (i=0; i<regionName.length() && i<FixedStringLengthForStore-1; ++i)
937 temp[i]=regionName[(G4int)i]; << 905 temp[i]=regionName[i];
938 fOut.write(temp, FixedStringLengthForSto 906 fOut.write(temp, FixedStringLengthForStore);
939 907
940 // cut values 908 // cut values
941 for (std::size_t idx=0; idx< NumberOfG4C 909 for (std::size_t idx=0; idx< NumberOfG4CutIndex; ++idx)
942 { 910 {
943 fOut.write( (char*)(&(cutValues[idx]) 911 fOut.write( (char*)(&(cutValues[idx])), sizeof(G4double));
944 } 912 }
945 } 913 }
946 } 914 }
947 fOut.close(); 915 fOut.close();
948 return true; 916 return true;
949 } 917 }
950 918
951 // ------------------------------------------- 919 // --------------------------------------------------------------------
952 G4bool 920 G4bool
953 G4ProductionCutsTable::CheckMaterialCutsCouple 921 G4ProductionCutsTable::CheckMaterialCutsCoupleInfo(const G4String& directory,
954 922 G4bool ascii)
955 { 923 {
956 // Check stored materialCutsCouple is consis 924 // Check stored materialCutsCouple is consistent
957 // with the current detector setup. 925 // with the current detector setup.
958 926
959 const G4String fileName = directory + "/" + 927 const G4String fileName = directory + "/" + "couple.dat";
960 const G4String key = "COUPLE-V3.0"; 928 const G4String key = "COUPLE-V3.0";
961 std::ifstream fIn; 929 std::ifstream fIn;
962 930
963 // open input file 931 // open input file
964 if (!ascii ) fIn.open(fileName,std::ios::in| 932 if (!ascii ) fIn.open(fileName,std::ios::in|std::ios::binary);
965 else fIn.open(fileName,std::ios::in) 933 else fIn.open(fileName,std::ios::in);
966 934
967 // check if the file has been opened success 935 // check if the file has been opened successfully
968 if (!fIn) 936 if (!fIn)
969 { 937 {
970 #ifdef G4VERBOSE 938 #ifdef G4VERBOSE
971 if (verboseLevel >0) 939 if (verboseLevel >0)
972 { 940 {
973 G4cerr << "G4ProductionCutTable::CheckMa 941 G4cerr << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo() - ";
974 G4cerr << "Cannot open file!" << fileNam 942 G4cerr << "Cannot open file!" << fileName << G4endl;
975 } 943 }
976 #endif 944 #endif
977 G4Exception( "G4ProductionCutsTable::Check 945 G4Exception( "G4ProductionCutsTable::CheckMaterialCutsCoupleInfo()",
978 "ProcCuts102", JustWarning, " 946 "ProcCuts102", JustWarning, "Cannot open file!");
979 return false; 947 return false;
980 } 948 }
981 949
982 char temp[FixedStringLengthForStore]; 950 char temp[FixedStringLengthForStore];
983 951
984 // key word 952 // key word
985 G4String keyword; 953 G4String keyword;
986 if (ascii) 954 if (ascii)
987 { 955 {
988 fIn >> keyword; 956 fIn >> keyword;
989 } 957 }
990 else 958 else
991 { 959 {
992 fIn.read(temp, FixedStringLengthForStore); 960 fIn.read(temp, FixedStringLengthForStore);
993 keyword = (const char*)(temp); 961 keyword = (const char*)(temp);
994 } 962 }
995 if (key!=keyword) 963 if (key!=keyword)
996 { 964 {
997 #ifdef G4VERBOSE 965 #ifdef G4VERBOSE
998 if (verboseLevel >0) 966 if (verboseLevel >0)
999 { 967 {
1000 G4cerr << "G4ProductionCutTable::CheckM 968 G4cerr << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo() - ";
1001 G4cerr << "Key word in " << fileName << 969 G4cerr << "Key word in " << fileName << "= " << keyword ;
1002 G4cerr <<"( should be "<< key << ")" 970 G4cerr <<"( should be "<< key << ")" << G4endl;
1003 } 971 }
1004 #endif 972 #endif
1005 G4Exception( "G4ProductionCutsTable::Chec 973 G4Exception( "G4ProductionCutsTable::CheckMaterialCutsCoupleInfo()",
1006 "ProcCuts103", JustWarning, 974 "ProcCuts103", JustWarning, "Bad Data Format");
1007 fIn.close(); 975 fIn.close();
1008 return false; 976 return false;
1009 } 977 }
1010 978
1011 // numberOfCouples 979 // numberOfCouples
1012 G4int numberOfCouples; 980 G4int numberOfCouples;
1013 if (ascii) 981 if (ascii)
1014 { 982 {
1015 fIn >> numberOfCouples; 983 fIn >> numberOfCouples;
1016 } 984 }
1017 else 985 else
1018 { 986 {
1019 fIn.read( (char*)(&numberOfCouples), size 987 fIn.read( (char*)(&numberOfCouples), sizeof(G4int));
1020 } 988 }
1021 989
1022 // Reset MCCIndexConversionTable 990 // Reset MCCIndexConversionTable
1023 mccConversionTable.Reset(numberOfCouples); 991 mccConversionTable.Reset(numberOfCouples);
1024 992
1025 // Read in couple information 993 // Read in couple information
1026 for (G4int idx=0; idx<numberOfCouples; ++id 994 for (G4int idx=0; idx<numberOfCouples; ++idx)
1027 { 995 {
1028 // read in index 996 // read in index
1029 G4int index; 997 G4int index;
1030 if (ascii) 998 if (ascii)
1031 { 999 {
1032 fIn >> index; 1000 fIn >> index;
1033 } 1001 }
1034 else 1002 else
1035 { 1003 {
1036 fIn.read( (char*)(&index), sizeof(G4int 1004 fIn.read( (char*)(&index), sizeof(G4int));
1037 } 1005 }
1038 // read in index material name 1006 // read in index material name
1039 char mat_name[FixedStringLengthForStore]; 1007 char mat_name[FixedStringLengthForStore];
1040 if (ascii) 1008 if (ascii)
1041 { 1009 {
1042 fIn >> mat_name; 1010 fIn >> mat_name;
1043 } 1011 }
1044 else 1012 else
1045 { 1013 {
1046 fIn.read(mat_name, FixedStringLengthFor 1014 fIn.read(mat_name, FixedStringLengthForStore);
1047 } 1015 }
1048 // read in index and region name 1016 // read in index and region name
1049 char region_name[FixedStringLengthForStor 1017 char region_name[FixedStringLengthForStore];
1050 if (ascii) 1018 if (ascii)
1051 { 1019 {
1052 fIn >> region_name; 1020 fIn >> region_name;
1053 } 1021 }
1054 else 1022 else
1055 { 1023 {
1056 fIn.read(region_name, FixedStringLength 1024 fIn.read(region_name, FixedStringLengthForStore);
1057 } 1025 }
1058 // cut value 1026 // cut value
1059 G4double cutValues[NumberOfG4CutIndex]; 1027 G4double cutValues[NumberOfG4CutIndex];
1060 for (std::size_t i=0; i< NumberOfG4CutInd 1028 for (std::size_t i=0; i< NumberOfG4CutIndex; ++i)
1061 { 1029 {
1062 if (ascii) 1030 if (ascii)
1063 { 1031 {
1064 fIn >> cutValues[i]; 1032 fIn >> cutValues[i];
1065 cutValues[i] *= (mm); 1033 cutValues[i] *= (mm);
1066 } 1034 }
1067 else 1035 else
1068 { 1036 {
1069 fIn.read( (char*)(&(cutValues[i])), s 1037 fIn.read( (char*)(&(cutValues[i])), sizeof(G4double));
1070 } 1038 }
1071 } 1039 }
1072 1040
1073 // Loop over all couples 1041 // Loop over all couples
1074 G4bool fOK = false; 1042 G4bool fOK = false;
1075 G4MaterialCutsCouple* aCouple = nullptr; 1043 G4MaterialCutsCouple* aCouple = nullptr;
1076 for (auto cItr=coupleTable.cbegin(); cItr 1044 for (auto cItr=coupleTable.cbegin(); cItr!=coupleTable.cend(); ++cItr)
1077 { 1045 {
1078 aCouple = (*cItr); 1046 aCouple = (*cItr);
1079 // check material name 1047 // check material name
1080 if ( mat_name != aCouple->GetMaterial( 1048 if ( mat_name != aCouple->GetMaterial()->GetName() ) continue;
1081 // check cut values 1049 // check cut values
1082 G4ProductionCuts* aCut = aCouple->GetPr 1050 G4ProductionCuts* aCut = aCouple->GetProductionCuts();
1083 G4bool fRatio = true; 1051 G4bool fRatio = true;
1084 for (std::size_t j=0; j< NumberOfG4CutI 1052 for (std::size_t j=0; j< NumberOfG4CutIndex; ++j)
1085 { 1053 {
1086 // check ratio only if values are not 1054 // check ratio only if values are not the same
1087 if (cutValues[j] != aCut->GetProducti << 1055 if (cutValues[j] != aCut->GetProductionCut(j))
1088 { 1056 {
1089 G4double ratio = cutValues[j]/aCut << 1057 G4double ratio = cutValues[j]/aCut->GetProductionCut(j);
1090 fRatio = fRatio && (0.999<ratio) && 1058 fRatio = fRatio && (0.999<ratio) && (ratio<1.001) ;
1091 } 1059 }
1092 } 1060 }
1093 if (!fRatio) continue; 1061 if (!fRatio) continue;
1094 // MCC matched 1062 // MCC matched
1095 fOK = true; 1063 fOK = true;
1096 mccConversionTable.SetNewIndex(index, a 1064 mccConversionTable.SetNewIndex(index, aCouple->GetIndex());
1097 break; 1065 break;
1098 } 1066 }
1099 1067
1100 if (fOK) 1068 if (fOK)
1101 { 1069 {
1102 #ifdef G4VERBOSE 1070 #ifdef G4VERBOSE
1103 // debug information 1071 // debug information
1104 if (verboseLevel >1) 1072 if (verboseLevel >1)
1105 { 1073 {
1106 G4String regionname(region_name); 1074 G4String regionname(region_name);
1107 G4Region* fRegion = nullptr; 1075 G4Region* fRegion = nullptr;
1108 if ( regionname != "NONE" ) 1076 if ( regionname != "NONE" )
1109 { 1077 {
1110 fRegion = fG4RegionStore->GetRegion 1078 fRegion = fG4RegionStore->GetRegion(region_name);
1111 if (fRegion == nullptr) 1079 if (fRegion == nullptr)
1112 { 1080 {
1113 G4cout << "G4ProductionCutTable:: 1081 G4cout << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo() - ";
1114 G4cout << "Region " << regionname 1082 G4cout << "Region " << regionname << " is not found ";
1115 G4cout << index << ": in " << fil 1083 G4cout << index << ": in " << fileName << G4endl;
1116 } 1084 }
1117 } 1085 }
1118 if (((regionname == "NONE") && (aCou 1086 if (((regionname == "NONE") && (aCouple->IsUsed()))
1119 || ((fRegion!=nullptr) && !IsCoupleU 1087 || ((fRegion!=nullptr) && !IsCoupleUsedInTheRegion(aCouple, fRegion)))
1120 { 1088 {
1121 G4cout << "G4ProductionCutTable::Ch 1089 G4cout << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo()"
1122 << G4endl; 1090 << G4endl;
1123 G4cout << "A Couple is used differe 1091 G4cout << "A Couple is used different region in the current setup ";
1124 G4cout << index << ": in " << fileN 1092 G4cout << index << ": in " << fileName << G4endl;
1125 G4cout << " material: " << mat_name 1093 G4cout << " material: " << mat_name ;
1126 G4cout << " region: " << region_nam 1094 G4cout << " region: " << region_name << G4endl;
1127 for (std::size_t ii=0; ii< NumberOf 1095 for (std::size_t ii=0; ii< NumberOfG4CutIndex; ++ii)
1128 { 1096 {
1129 G4cout << "cut[" << ii << "]=" << 1097 G4cout << "cut[" << ii << "]=" << cutValues[ii]/mm;
1130 G4cout << " mm : "; 1098 G4cout << " mm : ";
1131 } 1099 }
1132 G4cout << G4endl; 1100 G4cout << G4endl;
1133 } 1101 }
1134 else if ( index != aCouple->GetIndex 1102 else if ( index != aCouple->GetIndex() )
1135 { 1103 {
1136 G4cout << "G4ProductionCutTable::Ch 1104 G4cout << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo() - ";
1137 G4cout << "Index of couples was mod 1105 G4cout << "Index of couples was modified "<< G4endl;
1138 G4cout << aCouple->GetIndex() << ": 1106 G4cout << aCouple->GetIndex() << ":"
1139 << aCouple->GetMaterial()->G 1107 << aCouple->GetMaterial()->GetName();
1140 G4cout <<" is defined as " ; 1108 G4cout <<" is defined as " ;
1141 G4cout << index << ":" << mat_name 1109 G4cout << index << ":" << mat_name << " in " << fileName << G4endl;
1142 } 1110 }
1143 else 1111 else
1144 { 1112 {
1145 G4cout << "G4ProductionCutTable::Ch 1113 G4cout << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo() - ";
1146 G4cout << index << ":" << mat_name 1114 G4cout << index << ":" << mat_name << " in " << fileName ;
1147 G4cout << " is consistent with curr 1115 G4cout << " is consistent with current setup" << G4endl;
1148 } 1116 }
1149 } 1117 }
1150 #endif 1118 #endif
1151 } 1119 }
1152 else 1120 else
1153 { 1121 {
1154 #ifdef G4VERBOSE 1122 #ifdef G4VERBOSE
1155 if (verboseLevel >0) 1123 if (verboseLevel >0)
1156 { 1124 {
1157 G4cout << "G4ProductionCutTable::Chec 1125 G4cout << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo()"
1158 << G4endl; 1126 << G4endl;
1159 G4cout << "Couples are not defined in 1127 G4cout << "Couples are not defined in the current detector setup ";
1160 G4cout << index << ": in " << fileNam 1128 G4cout << index << ": in " << fileName << G4endl;
1161 G4cout << " material: " << mat_name ; 1129 G4cout << " material: " << mat_name ;
1162 G4cout << " region: " << region_name 1130 G4cout << " region: " << region_name << G4endl;
1163 for (std::size_t ii=0; ii< NumberOfG4 1131 for (std::size_t ii=0; ii< NumberOfG4CutIndex; ++ii)
1164 { 1132 {
1165 G4cout << "cut[" << ii << "]=" << c 1133 G4cout << "cut[" << ii << "]=" << cutValues[ii]/mm;
1166 G4cout << " mm : "; 1134 G4cout << " mm : ";
1167 } 1135 }
1168 G4cout << G4endl; 1136 G4cout << G4endl;
1169 } 1137 }
1170 #endif 1138 #endif
1171 } 1139 }
1172 } 1140 }
1173 fIn.close(); 1141 fIn.close();
1174 return true; 1142 return true;
1175 } 1143 }
1176 1144
1177 // ------------------------------------------ 1145 // --------------------------------------------------------------------
1178 G4bool G4ProductionCutsTable::StoreCutsInfo(c 1146 G4bool G4ProductionCutsTable::StoreCutsInfo(const G4String& directory,
1179 G 1147 G4bool ascii)
1180 { 1148 {
1181 // Store cut values information in files un 1149 // Store cut values information in files under the specified directory
1182 1150
1183 const G4String fileName = directory + "/" + 1151 const G4String fileName = directory + "/" + "cut.dat";
1184 const G4String key = "CUT-V3.0"; 1152 const G4String key = "CUT-V3.0";
1185 std::ofstream fOut; 1153 std::ofstream fOut;
1186 char temp[FixedStringLengthForStore]; 1154 char temp[FixedStringLengthForStore];
1187 1155
1188 // open output file 1156 // open output file
1189 if (!ascii ) fOut.open(fileName,std::ios::o 1157 if (!ascii ) fOut.open(fileName,std::ios::out|std::ios::binary);
1190 else fOut.open(fileName,std::ios::o 1158 else fOut.open(fileName,std::ios::out);
1191 1159
1192 // check if the file has been opened succes 1160 // check if the file has been opened successfully
1193 if (!fOut) 1161 if (!fOut)
1194 { 1162 {
1195 if(verboseLevel>0) 1163 if(verboseLevel>0)
1196 { 1164 {
1197 G4cerr << "G4ProductionCutsTable::Store 1165 G4cerr << "G4ProductionCutsTable::StoreCutsInfo() - ";
1198 G4cerr << "Cannot open file: " << fileN 1166 G4cerr << "Cannot open file: " << fileName << G4endl;
1199 } 1167 }
1200 G4Exception( "G4ProductionCutsTable::Stor 1168 G4Exception( "G4ProductionCutsTable::StoreCutsInfo()",
1201 "ProcCuts102", JustWarning, 1169 "ProcCuts102", JustWarning, "Cannot open file!");
1202 return false; 1170 return false;
1203 } 1171 }
1204 1172
1205 G4int numberOfCouples = (G4int)coupleTable. << 1173 G4int numberOfCouples = coupleTable.size();
1206 if (ascii) 1174 if (ascii)
1207 { 1175 {
1208 /////////////// ASCII mode ///////////// 1176 /////////////// ASCII mode /////////////////
1209 // key word 1177 // key word
1210 fOut << key << G4endl; 1178 fOut << key << G4endl;
1211 1179
1212 // number of couples in the table 1180 // number of couples in the table
1213 fOut << numberOfCouples << G4endl; 1181 fOut << numberOfCouples << G4endl;
1214 } 1182 }
1215 else 1183 else
1216 { 1184 {
1217 /////////////// Binary mode //////////// 1185 /////////////// Binary mode /////////////////
1218 // key word 1186 // key word
1219 std::size_t i; 1187 std::size_t i;
1220 for (i=0; i<FixedStringLengthForStore; ++ 1188 for (i=0; i<FixedStringLengthForStore; ++i)
1221 temp[i] = '\0'; 1189 temp[i] = '\0';
1222 for (i=0; i<key.length() && i<FixedString 1190 for (i=0; i<key.length() && i<FixedStringLengthForStore-1; ++i)
1223 temp[i]=key[(G4int)i]; << 1191 temp[i]=key[i];
1224 fOut.write(temp, FixedStringLengthForStor 1192 fOut.write(temp, FixedStringLengthForStore);
1225 1193
1226 // number of couples in the table 1194 // number of couples in the table
1227 fOut.write( (char*)(&numberOfCouples), si 1195 fOut.write( (char*)(&numberOfCouples), sizeof(G4int));
1228 } 1196 }
1229 1197
1230 for (std::size_t idx=0; idx <NumberOfG4CutI 1198 for (std::size_t idx=0; idx <NumberOfG4CutIndex; ++idx)
1231 { 1199 {
1232 const std::vector<G4double>* fRange = Ge 1200 const std::vector<G4double>* fRange = GetRangeCutsVector(idx);
1233 const std::vector<G4double>* fEnergy = Ge 1201 const std::vector<G4double>* fEnergy = GetEnergyCutsVector(idx);
1234 std::size_t i=0; 1202 std::size_t i=0;
1235 // Loop over all couples 1203 // Loop over all couples
1236 for (auto cItr=coupleTable.cbegin();cItr! 1204 for (auto cItr=coupleTable.cbegin();cItr!=coupleTable.cend(); ++cItr, ++i)
1237 { 1205 {
1238 if (ascii) 1206 if (ascii)
1239 { 1207 {
1240 /////////////// ASCII mode ///////// 1208 /////////////// ASCII mode /////////////////
1241 fOut.setf(std::ios::scientific); 1209 fOut.setf(std::ios::scientific);
1242 fOut << std::setw(20) << (*fRange)[i] 1210 fOut << std::setw(20) << (*fRange)[i]/mm ;
1243 fOut << std::setw(20) << (*fEnergy)[i 1211 fOut << std::setw(20) << (*fEnergy)[i]/keV << G4endl;
1244 fOut.unsetf(std::ios::scientific); 1212 fOut.unsetf(std::ios::scientific);
1245 } 1213 }
1246 else 1214 else
1247 { 1215 {
1248 /////////////// Binary mode //////// 1216 /////////////// Binary mode /////////////////
1249 G4double cut = (*fRange)[i]; 1217 G4double cut = (*fRange)[i];
1250 fOut.write((char*)(&cut), sizeof(G4do 1218 fOut.write((char*)(&cut), sizeof(G4double));
1251 cut = (*fEnergy)[i]; 1219 cut = (*fEnergy)[i];
1252 fOut.write((char*)(&cut), sizeof(G4do 1220 fOut.write((char*)(&cut), sizeof(G4double));
1253 } 1221 }
1254 } 1222 }
1255 } 1223 }
1256 fOut.close(); 1224 fOut.close();
1257 return true; 1225 return true;
1258 } 1226 }
1259 1227
1260 // ------------------------------------------ 1228 // --------------------------------------------------------------------
1261 G4bool G4ProductionCutsTable::RetrieveCutsInf 1229 G4bool G4ProductionCutsTable::RetrieveCutsInfo(const G4String& directory,
1262 1230 G4bool ascii)
1263 { 1231 {
1264 // Retrieve cut values information in files 1232 // Retrieve cut values information in files under the specified directory
1265 1233
1266 const G4String fileName = directory + "/" + 1234 const G4String fileName = directory + "/" + "cut.dat";
1267 const G4String key = "CUT-V3.0"; 1235 const G4String key = "CUT-V3.0";
1268 std::ifstream fIn; 1236 std::ifstream fIn;
1269 1237
1270 // open input file 1238 // open input file
1271 if (!ascii ) fIn.open(fileName,std::ios::in 1239 if (!ascii ) fIn.open(fileName,std::ios::in|std::ios::binary);
1272 else fIn.open(fileName,std::ios::in 1240 else fIn.open(fileName,std::ios::in);
1273 1241
1274 // check if the file has been opened succes 1242 // check if the file has been opened successfully
1275 if (!fIn) 1243 if (!fIn)
1276 { 1244 {
1277 if (verboseLevel >0) 1245 if (verboseLevel >0)
1278 { 1246 {
1279 G4cerr << "G4ProductionCutTable::Retrie 1247 G4cerr << "G4ProductionCutTable::RetrieveCutsInfo() - ";
1280 G4cerr << "Cannot open file: " << fileN 1248 G4cerr << "Cannot open file: " << fileName << G4endl;
1281 } 1249 }
1282 G4Exception( "G4ProductionCutsTable::Retr 1250 G4Exception( "G4ProductionCutsTable::RetrieveCutsInfo()",
1283 "ProcCuts102", JustWarning, 1251 "ProcCuts102", JustWarning, "Cannot open file!");
1284 return false; 1252 return false;
1285 } 1253 }
1286 1254
1287 char temp[FixedStringLengthForStore]; 1255 char temp[FixedStringLengthForStore];
1288 1256
1289 // key word 1257 // key word
1290 G4String keyword; 1258 G4String keyword;
1291 if (ascii) 1259 if (ascii)
1292 { 1260 {
1293 fIn >> keyword; 1261 fIn >> keyword;
1294 } 1262 }
1295 else 1263 else
1296 { 1264 {
1297 fIn.read(temp, FixedStringLengthForStore) 1265 fIn.read(temp, FixedStringLengthForStore);
1298 keyword = (const char*)(temp); 1266 keyword = (const char*)(temp);
1299 } 1267 }
1300 if (key!=keyword) 1268 if (key!=keyword)
1301 { 1269 {
1302 if (verboseLevel >0) 1270 if (verboseLevel >0)
1303 { 1271 {
1304 G4cerr << "G4ProductionCutTable::Retrie 1272 G4cerr << "G4ProductionCutTable::RetrieveCutsInfo() - ";
1305 G4cerr << "Key word in " << fileName << 1273 G4cerr << "Key word in " << fileName << "= " << keyword ;
1306 G4cerr <<"( should be "<< key << ")" 1274 G4cerr <<"( should be "<< key << ")" << G4endl;
1307 } 1275 }
1308 G4Exception( "G4ProductionCutsTable::Retr 1276 G4Exception( "G4ProductionCutsTable::RetrieveCutsInfo()",
1309 "ProcCuts103", JustWarning, 1277 "ProcCuts103", JustWarning, "Bad Data Format");
1310 return false; 1278 return false;
1311 } 1279 }
1312 1280
1313 // numberOfCouples 1281 // numberOfCouples
1314 G4int numberOfCouples; 1282 G4int numberOfCouples;
1315 if (ascii) 1283 if (ascii)
1316 { 1284 {
1317 fIn >> numberOfCouples; 1285 fIn >> numberOfCouples;
1318 if (fIn.fail()) 1286 if (fIn.fail())
1319 { 1287 {
1320 G4Exception( "G4ProductionCutsTable::Re 1288 G4Exception( "G4ProductionCutsTable::RetrieveCutsInfo()",
1321 "ProcCuts103", JustWarning 1289 "ProcCuts103", JustWarning, "Bad Data Format");
1322 return false; 1290 return false;
1323 } 1291 }
1324 } 1292 }
1325 else 1293 else
1326 { 1294 {
1327 fIn.read( (char*)(&numberOfCouples), size 1295 fIn.read( (char*)(&numberOfCouples), sizeof(G4int));
1328 } 1296 }
1329 1297
1330 if (numberOfCouples > static_cast<G4int>(mc 1298 if (numberOfCouples > static_cast<G4int>(mccConversionTable.size()) )
1331 { 1299 {
1332 G4Exception( "G4ProductionCutsTable::Retr 1300 G4Exception( "G4ProductionCutsTable::RetrieveCutsInfo()",
1333 "ProcCuts109", JustWarning, 1301 "ProcCuts109", JustWarning,
1334 "Number of Couples in the fi 1302 "Number of Couples in the file exceeds defined couples");
1335 } 1303 }
1336 numberOfCouples = (G4int)mccConversionTable << 1304 numberOfCouples = mccConversionTable.size();
1337 1305
1338 for (std::size_t idx=0; static_cast<G4int>( 1306 for (std::size_t idx=0; static_cast<G4int>(idx) <NumberOfG4CutIndex; ++idx)
1339 { 1307 {
1340 std::vector<G4double>* fRange = rangeCut 1308 std::vector<G4double>* fRange = rangeCutTable[idx];
1341 std::vector<G4double>* fEnergy = energyCu 1309 std::vector<G4double>* fEnergy = energyCutTable[idx];
1342 fRange->clear(); 1310 fRange->clear();
1343 fEnergy->clear(); 1311 fEnergy->clear();
1344 1312
1345 // Loop over all couples 1313 // Loop over all couples
1346 for (std::size_t i=0; static_cast<G4int>( 1314 for (std::size_t i=0; static_cast<G4int>(i)< numberOfCouples; ++i)
1347 { 1315 {
1348 G4double rcut, ecut; 1316 G4double rcut, ecut;
1349 if (ascii) 1317 if (ascii)
1350 { 1318 {
1351 fIn >> rcut >> ecut; 1319 fIn >> rcut >> ecut;
1352 if (fIn.fail()) 1320 if (fIn.fail())
1353 { 1321 {
1354 G4Exception( "G4ProductionCutsTable 1322 G4Exception( "G4ProductionCutsTable::RetrieveCutsInfo()",
1355 "ProcCuts103", JustWar 1323 "ProcCuts103", JustWarning, "Bad Data Format");
1356 return false; 1324 return false;
1357 } 1325 }
1358 rcut *= mm; 1326 rcut *= mm;
1359 ecut *= keV; 1327 ecut *= keV;
1360 } 1328 }
1361 else 1329 else
1362 { 1330 {
1363 fIn.read((char*)(&rcut), sizeof(G4dou 1331 fIn.read((char*)(&rcut), sizeof(G4double));
1364 fIn.read((char*)(&ecut), sizeof(G4dou 1332 fIn.read((char*)(&ecut), sizeof(G4double));
1365 } 1333 }
1366 if (!mccConversionTable.IsUsed(i)) cont 1334 if (!mccConversionTable.IsUsed(i)) continue;
1367 std::size_t new_index = mccConversionTa 1335 std::size_t new_index = mccConversionTable.GetIndex(i);
1368 (*fRange)[new_index] = rcut; 1336 (*fRange)[new_index] = rcut;
1369 (*fEnergy)[new_index] = ecut; 1337 (*fEnergy)[new_index] = ecut;
1370 } 1338 }
1371 } 1339 }
1372 return true; 1340 return true;
1373 } 1341 }
1374 1342
1375 // ------------------------------------------ 1343 // --------------------------------------------------------------------
1376 void G4ProductionCutsTable::SetVerboseLevel(G 1344 void G4ProductionCutsTable::SetVerboseLevel(G4int value)
1377 { 1345 {
1378 // Set same verbosity to all registered Ran 1346 // Set same verbosity to all registered RangeToEnergyConverters
1379 1347
1380 verboseLevel = value; 1348 verboseLevel = value;
1381 for (G4int ip=0; ip< NumberOfG4CutIndex; ++ 1349 for (G4int ip=0; ip< NumberOfG4CutIndex; ++ip)
1382 { 1350 {
1383 if (converters[ip] != nullptr ) 1351 if (converters[ip] != nullptr )
1384 { 1352 {
1385 converters[ip]->SetVerboseLevel(value); 1353 converters[ip]->SetVerboseLevel(value);
1386 } 1354 }
1387 } 1355 }
1388 } 1356 }
1389 1357
1390 // ------------------------------------------ 1358 // --------------------------------------------------------------------
1391 G4double G4ProductionCutsTable::GetMaxEnergyC 1359 G4double G4ProductionCutsTable::GetMaxEnergyCut()
1392 { 1360 {
1393 return G4VRangeToEnergyConverter::GetMaxEne 1361 return G4VRangeToEnergyConverter::GetMaxEnergyCut();
1394 } 1362 }
1395 1363
1396 // ------------------------------------------ 1364 // --------------------------------------------------------------------
1397 void G4ProductionCutsTable::SetMaxEnergyCut(G 1365 void G4ProductionCutsTable::SetMaxEnergyCut(G4double value)
1398 { 1366 {
1399 G4VRangeToEnergyConverter::SetMaxEnergyCut( 1367 G4VRangeToEnergyConverter::SetMaxEnergyCut(value);
1400 } 1368 }
1401 1369