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