61 if(fProductionCutsTable ==
nullptr)
63 fProductionCutsTable = &theProductionCutsTable;
65 return fProductionCutsTable;
73 rangeCutTable.push_back(
new std::vector<G4double>);
74 energyCutTable.push_back(
new std::vector<G4double>);
75 rangeDoubleVector[i] =
nullptr;
76 energyDoubleVector[i] =
nullptr;
77 converters[i] =
nullptr;
89 delete defaultProductionCuts;
90 defaultProductionCuts =
nullptr;
92 for(
auto itr=coupleTable.cbegin(); itr!=coupleTable.cend(); ++itr)
100 delete rangeCutTable[i];
101 delete energyCutTable[i];
102 delete converters[i];
103 if(rangeDoubleVector[i] !=
nullptr)
delete [] rangeDoubleVector[i];
104 if(energyDoubleVector[i] !=
nullptr)
delete [] energyDoubleVector[i];
105 rangeCutTable[i] =
nullptr;
106 energyCutTable[i] =
nullptr;
107 converters[i] =
nullptr;
108 rangeDoubleVector[i] =
nullptr;
109 energyDoubleVector[i] =
nullptr;
111 fProductionCutsTable =
nullptr;
114 fMessenger =
nullptr;
146 for(
auto CoupleItr=coupleTable.cbegin();
147 CoupleItr!=coupleTable.cend(); ++CoupleItr)
149 (*CoupleItr)->SetUseFlag(
false);
153 for(
auto rItr=fG4RegionStore->cbegin(); rItr!=fG4RegionStore->cend(); ++rItr)
159 if( (*rItr)->IsInMassGeometry() || (*rItr)->IsInParallelGeometry() )
162 auto mItr = (*rItr)->GetMaterialIterator();
163 std::size_t nMaterial = (*rItr)->GetNumberOfMaterials();
166 for(std::size_t iMate=0; iMate<nMaterial; ++iMate)
169 G4bool coupleAlreadyDefined =
false;
171 for(
auto cItr=coupleTable.cbegin(); cItr!=coupleTable.cend(); ++cItr)
173 if( (*cItr)->GetMaterial()==(*mItr)
174 && (*cItr)->GetProductionCuts()==fProductionCut)
176 coupleAlreadyDefined =
true;
183 if(!coupleAlreadyDefined)
186 coupleTable.push_back(aCouple);
187 aCouple->
SetIndex(coupleTable.size()-1);
191 (*rItr)->RegisterMaterialCouplePair((*mItr),aCouple);
196 auto rootLVItr = (*rItr)->GetRootLogicalVolumeIterator();
197 std::size_t nRootLV = (*rItr)->GetNumberOfRootVolumes();
198 for(std::size_t iLV=0; iLV<nRootLV; ++iLV)
204 ScanAndSetCouple(aLV,aCouple,aR);
220 std::size_t nCouple = coupleTable.size();
221 std::size_t nTable = energyCutTable[0]->size();
222 G4bool newCoupleAppears = nCouple>nTable;
225 for(std::size_t n=nCouple-nTable; n>0; --n)
229 rangeCutTable[nn]->push_back(-1.);
230 energyCutTable[nn]->push_back(-1.);
242 for(
auto cItr=coupleTable.cbegin(); cItr!=coupleTable.cend(); ++cItr)
246 if((*cItr)->IsRecalcNeeded())
251 (*(rangeCutTable[ptcl]))[idx] = rCut;
255 (*(energyCutTable[ptcl]))[idx] = converters[ptcl]->
Convert(rCut,aMat);
259 (*(energyCutTable[ptcl]))[idx] = -1.;
268 G4cout <<
"G4ProductionCutsTable::UpdateCoupleTable() - "
269 <<
"Elapsed time for calculation of energy cuts: " <<
G4endl;
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()];
290 for(std::size_t ixx=0; ixx<(*(rangeCutTable[ix])).size(); ++ixx)
292 rangeDoubleVector[ix][ixx] = (*(rangeCutTable[ix]))[ixx];
293 energyDoubleVector[ix][ixx] = (*(energyCutTable[ix]))[ixx];
313 ed <<
"Invoked prematurely before it is fully initialized.";
314 G4Exception(
"G4ProductionCutsTable::ConvertRangeToEnergy()",
322 if (material ==
nullptr)
return -1.0;
325 if (range == 0.0)
return 0.0;
326 if (range <0.0)
return -1.0;
331 if (index<0 || converters[index] ==
nullptr)
338 if(particle !=
nullptr)
341 { ed <<
"without valid particle pointer."; }
342 G4Exception(
"G4ProductionCutsTable::ConvertRangeToEnergy()",
349 return converters[index]->
Convert(range, material);
357 if (converters[i]!=0) converters[i]->
Reset();
385 if((aRegion!=
nullptr) && aLV->
GetRegion()!=aRegion)
return;
394 if(noDaughters==0)
return;
397 for(std::size_t i=0; i<noDaughters; ++i)
400 ScanAndSetCouple(daughterLVol,aCouple,aRegion);
408 G4cout <<
"========= Table of registered couples ============================"
410 for(
auto cItr=coupleTable.cbegin(); cItr!=coupleTable.cend(); ++cItr)
416 <<
" used in the geometry : ";
424 G4cout <<
" Range cuts : "
430 G4cout <<
" Energy thresholds : " ;
443 G4cout <<
" Region(s) which use this couple : " <<
G4endl;
444 for(
auto rItr=fG4RegionStore->cbegin();
445 rItr!=fG4RegionStore->cend(); ++rItr)
447 if (IsCoupleUsedInTheRegion(aCouple, *rItr) )
455 G4cout <<
"==================================================================" <<
G4endl;
472 G4cout <<
"G4ProductionCutsTable::StoreCutsTable()" <<
G4endl;
473 G4cout <<
" Material/Cuts information have been successfully stored ";
476 G4cout <<
" in Ascii mode ";
480 G4cout <<
" in Binary mode ";
497 G4cout <<
"G4ProductionCutsTable::RetrieveCutsTable()" <<
G4endl;
498 G4cout <<
" Material/Cuts information have been successfully retrieved ";
501 G4cout <<
" in Ascii mode ";
505 G4cout <<
" in Binary mode ";
521 G4cerr <<
"G4ProductionCutsTable::CheckForRetrieveCutsTable()"<<
G4endl;
526 G4cerr <<
"G4ProductionCutsTable::CheckMaterialInfo passed !!"<<
G4endl;
531 G4cerr <<
"G4ProductionCutsTable::CheckMaterialCutsCoupleInfo passed !!"
543 const G4String fileName = directory +
"/" +
"material.dat";
544 const G4String key =
"MATERIAL-V3.0";
548 if (!ascii ) fOut.open(fileName,std::ios::out|std::ios::binary);
549 else fOut.open(fileName,std::ios::out);
557 G4cerr <<
"G4ProductionCutsTable::StoreMaterialInfo() - ";
561 G4Exception(
"G4ProductionCutsTable::StoreMaterialInfo()",
568 G4int numberOfMaterial = matTable->size();
577 fOut << numberOfMaterial <<
G4endl;
579 fOut.setf(std::ios::scientific);
582 for (std::size_t idx=0;
static_cast<G4int>(idx)<numberOfMaterial; ++idx)
584 fOut << std::setw(FixedStringLengthForStore)
585 << ((*matTable)[idx])->GetName();
586 fOut << std::setw(FixedStringLengthForStore)
587 << ((*matTable)[idx])->GetDensity()/(g/cm3) <<
G4endl;
590 fOut.unsetf(std::ios::scientific);
596 char temp[FixedStringLengthForStore];
600 for (i=0; i<FixedStringLengthForStore; ++i)
604 for (i=0; i<key.length() && i<FixedStringLengthForStore-1; ++i)
608 fOut.write(temp, FixedStringLengthForStore);
611 fOut.write( (
char*)(&numberOfMaterial),
sizeof(
G4int));
614 for (std::size_t imat=0;
static_cast<G4int>(imat)<numberOfMaterial; ++imat)
616 G4String name = ((*matTable)[imat])->GetName();
617 G4double density = ((*matTable)[imat])->GetDensity();
618 for (i=0; i<FixedStringLengthForStore; ++i)
620 for (i=0; i<name.length() && i<FixedStringLengthForStore-1; ++i)
622 fOut.write(temp, FixedStringLengthForStore);
623 fOut.write( (
char*)(&density),
sizeof(
G4double));
637 const G4String fileName = directory +
"/" +
"material.dat";
638 const G4String key =
"MATERIAL-V3.0";
642 if (!ascii ) fIn.open(fileName,std::ios::in|std::ios::binary);
643 else fIn.open(fileName,std::ios::in);
651 G4cerr <<
"G4ProductionCutsTable::CheckMaterialInfo() - ";
655 G4Exception(
"G4ProductionCutsTable::CheckMaterialInfo()",
660 char temp[FixedStringLengthForStore];
670 fIn.read(temp, FixedStringLengthForStore);
671 keyword = (
const char*)(temp);
678 G4cerr <<
"G4ProductionCutsTable::CheckMaterialInfo() - ";
679 G4cerr <<
"Key word in " << fileName <<
"= " << keyword ;
683 G4Exception(
"G4ProductionCutsTable::CheckMaterialInfo()",
696 fIn.read( (
char*)(&nmat),
sizeof(
G4int));
698 if ((nmat<=0) || (nmat >100000))
700 G4Exception(
"G4ProductionCutsTable::CheckMaterialInfo()",
702 "Number of materials is less than zero or too big");
707 for (
G4int idx=0; idx<nmat ; ++idx)
715 G4cout <<
"G4ProductionCutsTable::CheckMaterialInfo() - ";
716 G4cout <<
"Encountered End of File " ;
725 char name[FixedStringLengthForStore];
729 fIn >> name >> density;
735 fIn.read(name, FixedStringLengthForStore);
736 fIn.read((
char*)(&density),
sizeof(
G4double));
743 G4cerr <<
"G4ProductionCutsTable::CheckMaterialInfo() - ";
744 G4cerr <<
"Bad data format ";
748 G4Exception(
"G4ProductionCutsTable::CheckMaterialInfo()",
755 if (aMaterial ==
nullptr )
continue;
758 if ((0.999>ratio) || (ratio>1.001) )
763 G4cerr <<
"G4ProductionCutsTable::CheckMaterialInfo() - ";
767 G4cerr <<
"Density:" << std::setiosflags(std::ios::scientific)
768 << density / (g/cm3) ;
771 G4cerr << std::resetiosflags(std::ios::scientific);
774 G4Exception(
"G4ProductionCutsTable::CheckMaterialInfo()",
775 "ProcCuts104",
JustWarning,
"Inconsistent material density");
792 const G4String fileName = directory +
"/" +
"couple.dat";
795 char temp[FixedStringLengthForStore];
798 if (!ascii ) fOut.open(fileName,std::ios::out|std::ios::binary);
799 else fOut.open(fileName,std::ios::out);
808 G4cerr <<
"G4ProductionCutsTable::StoreMaterialCutsCoupleInfo() - ";
812 G4Exception(
"G4ProductionCutsTable::StoreMaterialCutsCoupleInfo()",
817 G4int numberOfCouples = coupleTable.size();
822 fOut << std::setw(FixedStringLengthForStore) << key <<
G4endl;
825 fOut << numberOfCouples <<
G4endl;
832 for (i=0; i<FixedStringLengthForStore; ++i)
834 for (i=0; i<key.length() && i<FixedStringLengthForStore-1; ++i)
836 fOut.write(temp, FixedStringLengthForStore);
839 fOut.write( (
char*)(&numberOfCouples),
sizeof(
G4int));
843 for (
auto cItr=coupleTable.cbegin(); cItr!=coupleTable.cend(); ++cItr)
859 for(
auto rItr=fG4RegionStore->cbegin();
860 rItr!=fG4RegionStore->cend(); ++rItr)
862 if (IsCoupleUsedInTheRegion(aCouple, *rItr))
864 regionName = (*rItr)->GetName();
877 fOut << std::setw(FixedStringLengthForStore) << materialName<<
G4endl;
880 fOut << std::setw(FixedStringLengthForStore) << regionName<<
G4endl;
882 fOut.setf(std::ios::scientific);
886 fOut << std::setw(FixedStringLengthForStore) << cutValues[idx]/(mm)
889 fOut.unsetf(std::ios::scientific);
896 fOut.write( (
char*)(&index),
sizeof(
G4int));
900 for (i=0; i<FixedStringLengthForStore; ++i)
902 for (i=0; i<materialName.length() && i<FixedStringLengthForStore-1; ++i)
903 temp[i]=materialName[i];
904 fOut.write(temp, FixedStringLengthForStore);
907 for (i=0; i<FixedStringLengthForStore; ++i)
909 for (i=0; i<regionName.length() && i<FixedStringLengthForStore-1; ++i)
910 temp[i]=regionName[i];
911 fOut.write(temp, FixedStringLengthForStore);
916 fOut.write( (
char*)(&(cutValues[idx])),
sizeof(
G4double));
932 const G4String fileName = directory +
"/" +
"couple.dat";
937 if (!ascii ) fIn.open(fileName,std::ios::in|std::ios::binary);
938 else fIn.open(fileName,std::ios::in);
946 G4cerr <<
"G4ProductionCutTable::CheckMaterialCutsCoupleInfo() - ";
950 G4Exception(
"G4ProductionCutsTable::CheckMaterialCutsCoupleInfo()",
955 char temp[FixedStringLengthForStore];
965 fIn.read(temp, FixedStringLengthForStore);
966 keyword = (
const char*)(temp);
973 G4cerr <<
"G4ProductionCutTable::CheckMaterialCutsCoupleInfo() - ";
974 G4cerr <<
"Key word in " << fileName <<
"= " << keyword ;
978 G4Exception(
"G4ProductionCutsTable::CheckMaterialCutsCoupleInfo()",
985 G4int numberOfCouples;
988 fIn >> numberOfCouples;
992 fIn.read( (
char*)(&numberOfCouples),
sizeof(
G4int));
996 mccConversionTable.
Reset(numberOfCouples);
999 for (
G4int idx=0; idx<numberOfCouples; ++idx)
1009 fIn.read( (
char*)(&index),
sizeof(
G4int));
1012 char mat_name[FixedStringLengthForStore];
1019 fIn.read(mat_name, FixedStringLengthForStore);
1022 char region_name[FixedStringLengthForStore];
1029 fIn.read(region_name, FixedStringLengthForStore);
1037 fIn >> cutValues[i];
1038 cutValues[i] *= (mm);
1042 fIn.read( (
char*)(&(cutValues[i])),
sizeof(
G4double));
1049 for (
auto cItr=coupleTable.cbegin(); cItr!=coupleTable.cend(); ++cItr)
1063 fRatio = fRatio && (0.999<ratio) && (ratio<1.001) ;
1066 if (!fRatio)
continue;
1077 if (verboseLevel >1)
1081 if ( regionname !=
"NONE" )
1083 fRegion = fG4RegionStore->
GetRegion(region_name);
1084 if (fRegion ==
nullptr)
1086 G4cout <<
"G4ProductionCutTable::CheckMaterialCutsCoupleInfo() - ";
1087 G4cout <<
"Region " << regionname <<
" is not found ";
1091 if (((regionname ==
"NONE") && (aCouple->
IsUsed()))
1092 || ((fRegion!=
nullptr) && !IsCoupleUsedInTheRegion(aCouple, fRegion)))
1094 G4cout <<
"G4ProductionCutTable::CheckMaterialCutsCoupleInfo()"
1096 G4cout <<
"A Couple is used different region in the current setup ";
1098 G4cout <<
" material: " << mat_name ;
1102 G4cout <<
"cut[" << ii <<
"]=" << cutValues[ii]/mm;
1107 else if ( index != aCouple->
GetIndex() )
1109 G4cout <<
"G4ProductionCutTable::CheckMaterialCutsCoupleInfo() - ";
1113 G4cout <<
" is defined as " ;
1114 G4cout << index <<
":" << mat_name <<
" in " << fileName <<
G4endl;
1118 G4cout <<
"G4ProductionCutTable::CheckMaterialCutsCoupleInfo() - ";
1119 G4cout << index <<
":" << mat_name <<
" in " << fileName ;
1120 G4cout <<
" is consistent with current setup" <<
G4endl;
1128 if (verboseLevel >0)
1130 G4cout <<
"G4ProductionCutTable::CheckMaterialCutsCoupleInfo()"
1132 G4cout <<
"Couples are not defined in the current detector setup ";
1134 G4cout <<
" material: " << mat_name ;
1138 G4cout <<
"cut[" << ii <<
"]=" << cutValues[ii]/mm;
1156 const G4String fileName = directory +
"/" +
"cut.dat";
1159 char temp[FixedStringLengthForStore];
1162 if (!ascii ) fOut.open(fileName,std::ios::out|std::ios::binary);
1163 else fOut.open(fileName,std::ios::out);
1170 G4cerr <<
"G4ProductionCutsTable::StoreCutsInfo() - ";
1173 G4Exception(
"G4ProductionCutsTable::StoreCutsInfo()",
1178 G4int numberOfCouples = coupleTable.size();
1186 fOut << numberOfCouples <<
G4endl;
1193 for (i=0; i<FixedStringLengthForStore; ++i)
1195 for (i=0; i<key.length() && i<FixedStringLengthForStore-1; ++i)
1197 fOut.write(temp, FixedStringLengthForStore);
1200 fOut.write( (
char*)(&numberOfCouples),
sizeof(
G4int));
1209 for (
auto cItr=coupleTable.cbegin();cItr!=coupleTable.cend(); ++cItr, ++i)
1214 fOut.setf(std::ios::scientific);
1215 fOut << std::setw(20) << (*fRange)[i]/mm ;
1216 fOut << std::setw(20) << (*fEnergy)[i]/keV <<
G4endl;
1217 fOut.unsetf(std::ios::scientific);
1223 fOut.write((
char*)(&cut),
sizeof(
G4double));
1224 cut = (*fEnergy)[i];
1225 fOut.write((
char*)(&cut),
sizeof(
G4double));
1239 const G4String fileName = directory +
"/" +
"cut.dat";
1244 if (!ascii ) fIn.open(fileName,std::ios::in|std::ios::binary);
1245 else fIn.open(fileName,std::ios::in);
1250 if (verboseLevel >0)
1252 G4cerr <<
"G4ProductionCutTable::RetrieveCutsInfo() - ";
1255 G4Exception(
"G4ProductionCutsTable::RetrieveCutsInfo()",
1260 char temp[FixedStringLengthForStore];
1270 fIn.read(temp, FixedStringLengthForStore);
1271 keyword = (
const char*)(temp);
1275 if (verboseLevel >0)
1277 G4cerr <<
"G4ProductionCutTable::RetrieveCutsInfo() - ";
1278 G4cerr <<
"Key word in " << fileName <<
"= " << keyword ;
1281 G4Exception(
"G4ProductionCutsTable::RetrieveCutsInfo()",
1287 G4int numberOfCouples;
1290 fIn >> numberOfCouples;
1293 G4Exception(
"G4ProductionCutsTable::RetrieveCutsInfo()",
1300 fIn.read( (
char*)(&numberOfCouples),
sizeof(
G4int));
1303 if (numberOfCouples >
static_cast<G4int>(mccConversionTable.
size()) )
1305 G4Exception(
"G4ProductionCutsTable::RetrieveCutsInfo()",
1307 "Number of Couples in the file exceeds defined couples");
1309 numberOfCouples = mccConversionTable.
size();
1313 std::vector<G4double>* fRange = rangeCutTable[idx];
1314 std::vector<G4double>* fEnergy = energyCutTable[idx];
1319 for (std::size_t i=0;
static_cast<G4int>(i)< numberOfCouples; ++i)
1324 fIn >> rcut >> ecut;
1327 G4Exception(
"G4ProductionCutsTable::RetrieveCutsInfo()",
1336 fIn.read((
char*)(&rcut),
sizeof(
G4double));
1337 fIn.read((
char*)(&ecut),
sizeof(
G4double));
1339 if (!mccConversionTable.
IsUsed(i))
continue;
1340 std::size_t new_index = mccConversionTable.
GetIndex(i);
1341 (*fRange)[new_index] = rcut;
1342 (*fEnergy)[new_index] = ecut;
1353 verboseLevel = value;
1356 if (converters[ip] !=
nullptr )
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::ostringstream G4ExceptionDescription
std::vector< G4Material * > G4MaterialTable
G4GLOB_DLL std::ostream G4cerr
G4GLOB_DLL std::ostream G4cout
size_t GetNoDaughters() const
G4Region * GetRegion() const
G4Material * GetMaterial() const
G4VPhysicalVolume * GetDaughter(const G4int i) const
void SetMaterialCutsCouple(G4MaterialCutsCouple *cuts)
void SetNewIndex(std::size_t index, std::size_t new_value)
G4bool IsUsed(std::size_t index) const
G4int GetIndex(std::size_t index) const
void Reset(std::size_t size)
const G4Material * GetMaterial() const
void SetUseFlag(G4bool flg=true)
G4ProductionCuts * GetProductionCuts() const
G4double GetDensity() const
static G4MaterialTable * GetMaterialTable()
const G4String & GetName() const
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
const G4String & GetParticleName() const
static G4ParticleTable * GetParticleTable()
const std::vector< G4double > * GetRangeCutsVector(std::size_t pcIdx) const
virtual G4bool RetrieveCutsInfo(const G4String &directory, G4bool ascii=false)
G4bool RetrieveCutsTable(const G4String &directory, G4bool ascii=false)
virtual G4bool StoreCutsInfo(const G4String &directory, G4bool ascii=false)
G4double GetLowEdgeEnergy() const
void SetMaxEnergyCut(G4double value)
void UpdateCoupleTable(G4VPhysicalVolume *currentWorld)
G4double GetMaxEnergyCut()
void SetVerboseLevel(G4int value)
virtual G4bool CheckMaterialInfo(const G4String &directory, G4bool ascii=false)
G4int GetVerboseLevel() const
virtual G4bool StoreMaterialInfo(const G4String &directory, G4bool ascii=false)
G4double GetHighEdgeEnergy() const
virtual ~G4ProductionCutsTable()
const std::vector< G4double > * GetEnergyCutsVector(std::size_t pcIdx) const
G4bool StoreCutsTable(const G4String &directory, G4bool ascii=false)
G4bool CheckForRetrieveCutsTable(const G4String &directory, G4bool ascii=false)
void SetEnergyRange(G4double lowedge, G4double highedge)
static G4ProductionCutsTable * GetProductionCutsTable()
virtual G4bool CheckMaterialCutsCoupleInfo(const G4String &directory, G4bool ascii=false)
G4double ConvertRangeToEnergy(const G4ParticleDefinition *particle, const G4Material *material, G4double range)
virtual G4bool StoreMaterialCutsCoupleInfo(const G4String &directory, G4bool ascii=false)
static G4int GetIndex(const G4String &name)
G4double GetProductionCut(G4int index) const
const std::vector< G4double > & GetProductionCuts() const
static G4RegionStore * GetInstance()
G4Region * GetRegion(const G4String &name, G4bool verbose=true) const
G4LogicalVolume * GetLogicalVolume() const
void SetVerboseLevel(G4int value)
virtual G4double Convert(G4double rangeCut, const G4Material *material)
static G4double GetMaxEnergyCut()
static void SetMaxEnergyCut(G4double value)
static G4double GetLowEdgeEnergy()
static void SetEnergyRange(G4double lowedge, G4double highedge)
static G4double GetHighEdgeEnergy()