00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #include "G4ProductionCutsTable.hh"
00037 #include "G4ProductionCuts.hh"
00038 #include "G4MCCIndexConversionTable.hh"
00039 #include "G4ProductionCutsTableMessenger.hh"
00040 #include "G4ParticleDefinition.hh"
00041 #include "G4ParticleTable.hh"
00042 #include "G4RegionStore.hh"
00043 #include "G4LogicalVolume.hh"
00044 #include "G4VPhysicalVolume.hh"
00045 #include "G4RToEConvForElectron.hh"
00046 #include "G4RToEConvForGamma.hh"
00047 #include "G4RToEConvForPositron.hh"
00048 #include "G4RToEConvForProton.hh"
00049 #include "G4MaterialTable.hh"
00050 #include "G4Material.hh"
00051 #include "G4UnitsTable.hh"
00052
00053 #include "G4Timer.hh"
00054 #include "G4SystemOfUnits.hh"
00055 #include "G4ios.hh"
00056 #include <iomanip>
00057 #include <fstream>
00058
00059
00060 G4ProductionCutsTable* G4ProductionCutsTable::fG4ProductionCutsTable = 0;
00061
00063 G4ProductionCutsTable* G4ProductionCutsTable::GetProductionCutsTable()
00064 {
00065 static G4ProductionCutsTable theProductionCutsTable;
00066 if(!fG4ProductionCutsTable){
00067 fG4ProductionCutsTable = &theProductionCutsTable;
00068 }
00069 return fG4ProductionCutsTable;
00070 }
00071
00073 G4ProductionCutsTable::G4ProductionCutsTable()
00074 : firstUse(true),verboseLevel(1),fMessenger(0)
00075 {
00076 for(size_t i=0;i< NumberOfG4CutIndex;i++)
00077 {
00078 rangeCutTable.push_back(new G4CutVectorForAParticle);
00079 energyCutTable.push_back(new G4CutVectorForAParticle);
00080 rangeDoubleVector[i] = 0;
00081 energyDoubleVector[i] = 0;
00082 converters[i] = 0;
00083 }
00084 fG4RegionStore = G4RegionStore::GetInstance();
00085 defaultProductionCuts = new G4ProductionCuts();
00086
00087
00088 fMessenger = new G4ProductionCutsTableMessenger(this);
00089 }
00090
00092 G4ProductionCutsTable::G4ProductionCutsTable(const G4ProductionCutsTable& )
00093 {
00094 fMessenger=0;
00095 defaultProductionCuts=0;
00096 fG4RegionStore =0;
00097 }
00098
00100 G4ProductionCutsTable::~G4ProductionCutsTable()
00101 {
00102 if (defaultProductionCuts !=0) {
00103 delete defaultProductionCuts;
00104 defaultProductionCuts =0;
00105 }
00106
00107 for(CoupleTableIterator itr=coupleTable.begin();itr!=coupleTable.end();itr++){
00108 delete (*itr);
00109 }
00110 coupleTable.clear();
00111
00112 for(size_t i=0;i< NumberOfG4CutIndex;i++){
00113 delete rangeCutTable[i];
00114 delete energyCutTable[i];
00115 delete converters[i];
00116 if(rangeDoubleVector[i]!=0) delete [] rangeDoubleVector[i];
00117 if(energyDoubleVector[i]!=0) delete [] energyDoubleVector[i];
00118 }
00119 fG4ProductionCutsTable =0;
00120
00121 if (fMessenger !=0) delete fMessenger;
00122 fMessenger = 0;
00123 }
00124
00125 void G4ProductionCutsTable::UpdateCoupleTable(G4VPhysicalVolume* )
00126 {
00127 if(firstUse){
00128 if(G4ParticleTable::GetParticleTable()->FindParticle("gamma")){
00129 converters[0] = new G4RToEConvForGamma();
00130 converters[0]->SetVerboseLevel(GetVerboseLevel());
00131 }
00132 if(G4ParticleTable::GetParticleTable()->FindParticle("e-")){
00133 converters[1] = new G4RToEConvForElectron();
00134 converters[1]->SetVerboseLevel(GetVerboseLevel());
00135 }
00136 if(G4ParticleTable::GetParticleTable()->FindParticle("e+")){
00137 converters[2] = new G4RToEConvForPositron();
00138 converters[2]->SetVerboseLevel(GetVerboseLevel());
00139 }
00140 if(G4ParticleTable::GetParticleTable()->FindParticle("proton")){
00141 converters[3] = new G4RToEConvForProton();
00142 converters[3]->SetVerboseLevel(GetVerboseLevel());
00143 }
00144 firstUse = false;
00145 }
00146
00147
00148 for(CoupleTableIterator CoupleItr=coupleTable.begin();
00149 CoupleItr!=coupleTable.end();CoupleItr++)
00150 {
00151 (*CoupleItr)->SetUseFlag(false);
00152 }
00153
00154
00155 typedef std::vector<G4Region*>::iterator regionIterator;
00156 for(regionIterator rItr=fG4RegionStore->begin();
00157 rItr!=fG4RegionStore->end();rItr++)
00158 {
00159
00160
00161
00162 if((*rItr)->IsInMassGeometry() || (*rItr)->IsInParallelGeometry())
00163 {
00164
00165 G4ProductionCuts* fProductionCut = (*rItr)->GetProductionCuts();
00166 std::vector<G4Material*>::const_iterator mItr =
00167 (*rItr)->GetMaterialIterator();
00168 size_t nMaterial = (*rItr)->GetNumberOfMaterials();
00169 (*rItr)->ClearMap();
00170
00171 for(size_t iMate=0;iMate<nMaterial;iMate++){
00172
00173 G4bool coupleAlreadyDefined = false;
00174 G4MaterialCutsCouple* aCouple;
00175 for(CoupleTableIterator cItr=coupleTable.begin();
00176 cItr!=coupleTable.end();cItr++){
00177 if( (*cItr)->GetMaterial()==(*mItr) &&
00178 (*cItr)->GetProductionCuts()==fProductionCut){
00179 coupleAlreadyDefined = true;
00180 aCouple = *cItr;
00181 break;
00182 }
00183 }
00184
00185
00186 if(!coupleAlreadyDefined){
00187 aCouple = new G4MaterialCutsCouple((*mItr),fProductionCut);
00188 coupleTable.push_back(aCouple);
00189 aCouple->SetIndex(coupleTable.size()-1);
00190 }
00191
00192
00193 (*rItr)->RegisterMaterialCouplePair((*mItr),aCouple);
00194
00195
00196 aCouple->SetUseFlag();
00197
00198 std::vector<G4LogicalVolume*>::iterator rootLVItr
00199 = (*rItr)->GetRootLogicalVolumeIterator();
00200 size_t nRootLV = (*rItr)->GetNumberOfRootVolumes();
00201 for(size_t iLV=0;iLV<nRootLV;iLV++) {
00202
00203 G4LogicalVolume* aLV = *rootLVItr;
00204 G4Region* aR = *rItr;
00205
00206 ScanAndSetCouple(aLV,aCouple,aR);
00207
00208
00209 rootLVItr++;
00210 }
00211
00212
00213 mItr++;
00214 }
00215 }
00216 }
00217
00218
00219
00220
00221
00222 size_t nCouple = coupleTable.size();
00223 size_t nTable = energyCutTable[0]->size();
00224 G4bool newCoupleAppears = nCouple>nTable;
00225 if(newCoupleAppears) {
00226 for(size_t n=nCouple-nTable;n>0;n--) {
00227 for(size_t nn=0;nn< NumberOfG4CutIndex;nn++){
00228 rangeCutTable[nn]->push_back(-1.);
00229 energyCutTable[nn]->push_back(-1.);
00230 }
00231 }
00232 }
00233
00234
00235 size_t idx = 0;
00236 G4Timer timer;
00237 if (verboseLevel>2) {
00238 timer.Start();
00239 }
00240 for(CoupleTableIterator cItr=coupleTable.begin();
00241 cItr!=coupleTable.end();cItr++){
00242 G4ProductionCuts* aCut = (*cItr)->GetProductionCuts();
00243 const G4Material* aMat = (*cItr)->GetMaterial();
00244 if((*cItr)->IsRecalcNeeded()) {
00245 for(size_t ptcl=0;ptcl< NumberOfG4CutIndex;ptcl++){
00246 G4double rCut = aCut->GetProductionCut(ptcl);
00247 (*(rangeCutTable[ptcl]))[idx] = rCut;
00248
00249 if(converters[ptcl]) {
00250 (*(energyCutTable[ptcl]))[idx] = converters[ptcl]->Convert(rCut,aMat);
00251 } else {
00252 (*(energyCutTable[ptcl]))[idx] = -1.;
00253 }
00254 }
00255 }
00256 idx++;
00257 }
00258 if (verboseLevel>2) {
00259 timer.Stop();
00260 G4cout << "G4ProductionCutsTable::UpdateCoupleTable "
00261 << " elapsed time for calculation of energy cuts " << G4endl;
00262 G4cout << timer <<G4endl;
00263 }
00264
00265
00266 if(newCoupleAppears){
00267 for(size_t ix=0;ix<NumberOfG4CutIndex;ix++){
00268 G4double* rangeVOld = rangeDoubleVector[ix];
00269 G4double* energyVOld = energyDoubleVector[ix];
00270 if(rangeVOld) delete [] rangeVOld;
00271 if(energyVOld) delete [] energyVOld;
00272 rangeDoubleVector[ix] = new G4double[(*(rangeCutTable[ix])).size()];
00273 energyDoubleVector[ix] = new G4double[(*(energyCutTable[ix])).size()];
00274 }
00275 }
00276
00277
00278 for(size_t ix=0;ix<NumberOfG4CutIndex;ix++) {
00279 for(size_t ixx=0;ixx<(*(rangeCutTable[ix])).size();ixx++) {
00280 rangeDoubleVector[ix][ixx] = (*(rangeCutTable[ix]))[ixx];
00281 energyDoubleVector[ix][ixx] = (*(energyCutTable[ix]))[ixx];
00282 }
00283 }
00284 }
00285
00286
00288 G4double G4ProductionCutsTable::ConvertRangeToEnergy(
00289 const G4ParticleDefinition* particle,
00290 const G4Material* material,
00291 G4double range
00292 )
00293 {
00294
00295
00296 if (material ==0) return -1.0;
00297
00298
00299 if (range ==0.0) return 0.0;
00300 if (range <0.0) return -1.0;
00301
00302
00303 G4int index = G4ProductionCuts::GetIndex(particle);
00304
00305 if (index<0) {
00306 #ifdef G4VERBOSE
00307 if (verboseLevel >1) {
00308 G4cout << "G4ProductionCutsTable::ConvertRangeToEnergy" ;
00309 G4cout << particle->GetParticleName() << " has no cut value " << G4endl;
00310 }
00311 #endif
00312 return -1.0;
00313 }
00314
00315 return converters[index]->Convert(range, material);
00316
00317 }
00318
00320 void G4ProductionCutsTable::ResetConverters()
00321 {
00322 for(size_t i=0;i< NumberOfG4CutIndex;i++){
00323 if (converters[i]!=0) converters[i]->Reset();
00324 }
00325 }
00326
00327
00329 void G4ProductionCutsTable::SetEnergyRange(G4double lowedge, G4double highedge)
00330 {
00331 G4VRangeToEnergyConverter::SetEnergyRange(lowedge,highedge);
00332 }
00333
00335 G4double G4ProductionCutsTable::GetLowEdgeEnergy() const
00336 {
00337 return G4VRangeToEnergyConverter::GetLowEdgeEnergy();
00338 }
00339
00341 G4double G4ProductionCutsTable::GetHighEdgeEnergy() const
00342 {
00343 return G4VRangeToEnergyConverter::GetHighEdgeEnergy();
00344 }
00345
00346
00348 void G4ProductionCutsTable::ScanAndSetCouple(G4LogicalVolume* aLV,
00349 G4MaterialCutsCouple* aCouple,
00350 G4Region* aRegion)
00351 {
00352
00353 if((aRegion!=0) && aLV->GetRegion()!=aRegion) return;
00354
00355
00356 if(aLV->GetMaterial()==aCouple->GetMaterial()) {
00357 aLV->SetMaterialCutsCouple(aCouple);
00358 }
00359
00360 size_t noDaughters = aLV->GetNoDaughters();
00361 if(noDaughters==0) return;
00362
00363
00364 for(size_t i=0;i<noDaughters;i++){
00365 G4LogicalVolume* daughterLVol = aLV->GetDaughter(i)->GetLogicalVolume();
00366 ScanAndSetCouple(daughterLVol,aCouple,aRegion);
00367 }
00368 }
00369
00371 void G4ProductionCutsTable::DumpCouples() const
00372 {
00373 G4cout << G4endl;
00374 G4cout << "========= Table of registered couples =============================="
00375 << G4endl;
00376 for(CoupleTableIterator cItr=coupleTable.begin();
00377 cItr!=coupleTable.end();cItr++) {
00378 G4MaterialCutsCouple* aCouple = (*cItr);
00379 G4ProductionCuts* aCut = aCouple->GetProductionCuts();
00380 G4cout << G4endl;
00381 G4cout << "Index : " << aCouple->GetIndex()
00382 << " used in the geometry : ";
00383 if(aCouple->IsUsed()) G4cout << "Yes";
00384 else G4cout << "No ";
00385 G4cout << " recalculation needed : ";
00386 if(aCouple->IsRecalcNeeded()) G4cout << "Yes";
00387 else G4cout << "No ";
00388 G4cout << G4endl;
00389 G4cout << " Material : " << aCouple->GetMaterial()->GetName() << G4endl;
00390 G4cout << " Range cuts : "
00391 << " gamma " << G4BestUnit(aCut->GetProductionCut("gamma"),"Length")
00392 << " e- " << G4BestUnit(aCut->GetProductionCut("e-"),"Length")
00393 << " e+ " << G4BestUnit(aCut->GetProductionCut("e+"),"Length")
00394 << " proton " << G4BestUnit(aCut->GetProductionCut("proton"),"Length")
00395 << G4endl;
00396 G4cout << " Energy thresholds : " ;
00397 if(aCouple->IsRecalcNeeded()) {
00398 G4cout << " is not ready to print";
00399 } else {
00400 G4cout << " gamma " << G4BestUnit((*(energyCutTable[0]))[aCouple->GetIndex()],"Energy")
00401 << " e- " << G4BestUnit((*(energyCutTable[1]))[aCouple->GetIndex()],"Energy")
00402 << " e+ " << G4BestUnit((*(energyCutTable[2]))[aCouple->GetIndex()],"Energy")
00403 << " proton " << G4BestUnit((*(energyCutTable[3]))[aCouple->GetIndex()],"Energy");
00404 }
00405 G4cout << G4endl;
00406
00407 if(aCouple->IsUsed()) {
00408 G4cout << " Region(s) which use this couple : " << G4endl;
00409 typedef std::vector<G4Region*>::iterator regionIterator;
00410 for(regionIterator rItr=fG4RegionStore->begin();
00411 rItr!=fG4RegionStore->end();rItr++) {
00412 if (IsCoupleUsedInTheRegion(aCouple, *rItr) ){
00413 G4cout << " " << (*rItr)->GetName() << G4endl;
00414 }
00415 }
00416 }
00417 }
00418 G4cout << G4endl;
00419 G4cout << "====================================================================" << G4endl;
00420 G4cout << G4endl;
00421 }
00422
00423
00425
00426 G4bool G4ProductionCutsTable::StoreCutsTable(const G4String& dir,
00427 G4bool ascii)
00428 {
00429 if (!StoreMaterialInfo(dir, ascii)) return false;
00430 if (!StoreMaterialCutsCoupleInfo(dir, ascii)) return false;
00431 if (!StoreCutsInfo(dir, ascii)) return false;
00432
00433 #ifdef G4VERBOSE
00434 if (verboseLevel >2) {
00435 G4cout << "G4ProductionCutsTable::StoreCutsTable " ;
00436 G4cout << " Material/Cuts information have been succesfully stored ";
00437 if (ascii) {
00438 G4cout << " in Ascii mode ";
00439 }else{
00440 G4cout << " in Binary mode ";
00441 }
00442 G4cout << " under " << dir << G4endl;
00443 }
00444 #endif
00445 return true;
00446 }
00447
00449 G4bool G4ProductionCutsTable::RetrieveCutsTable(const G4String& dir,
00450 G4bool ascii)
00451 {
00452 if (!CheckForRetrieveCutsTable(dir, ascii)) return false;
00453 if (!RetrieveCutsInfo(dir, ascii)) return false;
00454 #ifdef G4VERBOSE
00455 if (verboseLevel >2) {
00456 G4cout << "G4ProductionCutsTable::RetrieveCutsTable " ;
00457 G4cout << " Material/Cuts information have been succesfully retreived ";
00458 if (ascii) {
00459 G4cout << " in Ascii mode ";
00460 }else{
00461 G4cout << " in Binary mode ";
00462 }
00463 G4cout << " under " << dir << G4endl;
00464 }
00465 #endif
00466 return true;
00467 }
00468
00470
00471
00472
00473 G4bool
00474 G4ProductionCutsTable::CheckForRetrieveCutsTable(const G4String& directory,
00475 G4bool ascii)
00476 {
00477 G4cerr << "G4ProductionCutsTable::CheckForRetrieveCutsTable!!"<< G4endl;
00478
00479 if (!CheckMaterialInfo(directory, ascii)) return false;
00480 if (verboseLevel >2) {
00481 G4cerr << "G4ProductionCutsTable::CheckMaterialInfo passed !!"<< G4endl;
00482 }
00483 if (!CheckMaterialCutsCoupleInfo(directory, ascii)) return false;
00484 if (verboseLevel >2) {
00485 G4cerr << "G4ProductionCutsTable::CheckMaterialCutsCoupleInfo passed !!"<< G4endl;
00486 }
00487 return true;
00488 }
00489
00491
00492
00493 G4bool G4ProductionCutsTable::StoreMaterialInfo(const G4String& directory,
00494 G4bool ascii)
00495 {
00496 const G4String fileName = directory + "/" + "material.dat";
00497 const G4String key = "MATERIAL-V3.0";
00498 std::ofstream fOut;
00499
00500
00501 if (!ascii )
00502 fOut.open(fileName,std::ios::out|std::ios::binary);
00503 else
00504 fOut.open(fileName,std::ios::out);
00505
00506
00507
00508 if (!fOut) {
00509 #ifdef G4VERBOSE
00510 if (verboseLevel>0) {
00511 G4cerr << "G4ProductionCutsTable::StoreMaterialInfo ";
00512 G4cerr << " Can not open file " << fileName << G4endl;
00513 }
00514 #endif
00515 G4Exception( "G4ProductionCutsTable::StoreMaterialInfo()",
00516 "ProcCuts102",
00517 JustWarning, "Can not open file ");
00518 return false;
00519 }
00520
00521 const G4MaterialTable* matTable = G4Material::GetMaterialTable();
00522
00523 G4int numberOfMaterial = matTable->size();
00524
00525 if (ascii) {
00527
00528 fOut << key << G4endl;
00529
00530
00531 fOut << numberOfMaterial << G4endl;
00532
00533 fOut.setf(std::ios::scientific);
00534
00535
00536 for (size_t idx=0; static_cast<G4int>(idx)<numberOfMaterial; ++idx){
00537 fOut << std::setw(FixedStringLengthForStore)
00538 << ((*matTable)[idx])->GetName();
00539 fOut << std::setw(FixedStringLengthForStore)
00540 << ((*matTable)[idx])->GetDensity()/(g/cm3) << G4endl;
00541 }
00542
00543 fOut.unsetf(std::ios::scientific);
00544
00545 } else {
00547 char temp[FixedStringLengthForStore];
00548 size_t i;
00549
00550
00551 for (i=0; i<FixedStringLengthForStore; ++i){
00552 temp[i] = '\0';
00553 }
00554 for (i=0; i<key.length() && i<FixedStringLengthForStore-1; ++i){
00555 temp[i]=key[i];
00556 }
00557 fOut.write(temp, FixedStringLengthForStore);
00558
00559
00560 fOut.write( (char*)(&numberOfMaterial), sizeof (G4int));
00561
00562
00563 for (size_t imat=0; static_cast<G4int>(imat)<numberOfMaterial; ++imat){
00564 G4String name = ((*matTable)[imat])->GetName();
00565 G4double density = ((*matTable)[imat])->GetDensity();
00566 for (i=0; i<FixedStringLengthForStore; ++i) temp[i] = '\0';
00567 for (i=0; i<name.length() && i<FixedStringLengthForStore-1; ++i)
00568 temp[i]=name[i];
00569 fOut.write(temp, FixedStringLengthForStore);
00570 fOut.write( (char*)(&density), sizeof (G4double));
00571 }
00572 }
00573
00574 fOut.close();
00575 return true;
00576 }
00577
00579
00580 G4bool G4ProductionCutsTable::CheckMaterialInfo(const G4String& directory,
00581 G4bool ascii)
00582 {
00583 const G4String fileName = directory + "/" + "material.dat";
00584 const G4String key = "MATERIAL-V3.0";
00585 std::ifstream fIn;
00586
00587
00588 if (!ascii )
00589 fIn.open(fileName,std::ios::in|std::ios::binary);
00590 else
00591 fIn.open(fileName,std::ios::in);
00592
00593
00594 if (!fIn) {
00595 #ifdef G4VERBOSE
00596 if (verboseLevel >0) {
00597 G4cerr << "G4ProductionCutsTable::CheckMaterialInfo ";
00598 G4cerr << " Can not open file " << fileName << G4endl;
00599 }
00600 #endif
00601 G4Exception( "G4ProductionCutsTable::CheckMaterialInfo()",
00602 "ProcCuts102",
00603 JustWarning, "Can not open file ");
00604 return false;
00605 }
00606
00607 char temp[FixedStringLengthForStore];
00608
00609
00610 G4String keyword;
00611 if (ascii) {
00612 fIn >> keyword;
00613 } else {
00614 fIn.read(temp, FixedStringLengthForStore);
00615 keyword = (const char*)(temp);
00616 }
00617 if (key!=keyword) {
00618 #ifdef G4VERBOSE
00619 if (verboseLevel >0) {
00620 G4cerr << "G4ProductionCutsTable::CheckMaterialInfo ";
00621 G4cerr << " Key word in " << fileName << "= " << keyword ;
00622 G4cerr <<"( should be "<< key << ")" <<G4endl;
00623 }
00624 #endif
00625 G4Exception( "G4ProductionCutsTable::CheckMaterialInfo()",
00626 "ProcCuts103",
00627 JustWarning, "Bad Data Format");
00628 return false;
00629 }
00630
00631
00632 G4int nmat;
00633 if (ascii) {
00634 fIn >> nmat;
00635 } else {
00636 fIn.read( (char*)(&nmat), sizeof (G4int));
00637 }
00638 if ((nmat<=0) || (nmat >100000)){
00639 G4Exception( "G4ProductionCutsTable::CheckMaterialInfo()",
00640 "ProcCuts108",JustWarning,
00641 "Number of materials is less than zero or too big");
00642 return false;
00643 }
00644
00645
00646 for (G4int idx=0; idx<nmat ; ++idx){
00647
00648 if(fIn.eof()) {
00649 #ifdef G4VERBOSE
00650 if (verboseLevel >0) {
00651 G4cout << "G4ProductionCutsTable::CheckMaterialInfo ";
00652 G4cout << " encountered End of File " ;
00653 G4cout << " at " << idx+1 << "th material "<< G4endl;
00654 }
00655 #endif
00656 fIn.close();
00657 return false;
00658 }
00659
00660
00661 char name[FixedStringLengthForStore];
00662 double density;
00663 if (ascii) {
00664 fIn >> name >> density;
00665 density *= (g/cm3);
00666
00667 } else {
00668 fIn.read(name, FixedStringLengthForStore);
00669 fIn.read((char*)(&density), sizeof (G4double));
00670 }
00671 if (fIn.fail()) {
00672 #ifdef G4VERBOSE
00673 if (verboseLevel >0) {
00674 G4cerr << "G4ProductionCutsTable::CheckMaterialInfo ";
00675 G4cerr << " Bad data format ";
00676 G4cerr << " at " << idx+1 << "th material "<< G4endl;
00677 }
00678 #endif
00679 G4Exception( "G4ProductionCutsTable::CheckMaterialInfo()",
00680 "ProcCuts103",
00681 JustWarning, "Bad Data Format");
00682 fIn.close();
00683 return false;
00684 }
00685
00686 G4Material* aMaterial = G4Material::GetMaterial(name);
00687 if (aMaterial ==0 ) continue;
00688
00689 G4double ratio = std::fabs(density/aMaterial->GetDensity() );
00690 if ((0.999>ratio) || (ratio>1.001) ){
00691 #ifdef G4VERBOSE
00692 if (verboseLevel >0) {
00693 G4cerr << "G4ProductionCutsTable::CheckMaterialInfo ";
00694 G4cerr << " Inconsistent material density" << G4endl;;
00695 G4cerr << " at " << idx+1 << "th material "<< G4endl;
00696 G4cerr << "Name: " << name << G4endl;
00697 G4cerr << "Density:" << std::setiosflags(std::ios::scientific) << density / (g/cm3) ;
00698 G4cerr << "(should be " << aMaterial->GetDensity()/(g/cm3)<< ")" << " [g/cm3]"<< G4endl;
00699 G4cerr << std::resetiosflags(std::ios::scientific);
00700 }
00701 #endif
00702 G4Exception( "G4ProductionCutsTable::CheckMaterialInfo()",
00703 "ProcCuts104",
00704 JustWarning, "Inconsitent matrial density");
00705 fIn.close();
00706 return false;
00707 }
00708
00709 }
00710
00711 fIn.close();
00712 return true;
00713
00714 }
00715
00716
00718
00719
00720 G4bool
00721 G4ProductionCutsTable::StoreMaterialCutsCoupleInfo(const G4String& directory,
00722 G4bool ascii)
00723 {
00724 const G4String fileName = directory + "/" + "couple.dat";
00725 const G4String key = "COUPLE-V3.0";
00726 std::ofstream fOut;
00727 char temp[FixedStringLengthForStore];
00728
00729
00730 if (!ascii )
00731 fOut.open(fileName,std::ios::out|std::ios::binary);
00732 else
00733 fOut.open(fileName,std::ios::out);
00734
00735
00736
00737 if (!fOut) {
00738 #ifdef G4VERBOSE
00739 if (verboseLevel >0) {
00740 G4cerr << "G4ProductionCutsTable::StoreMaterialCutsCoupleInfo ";
00741 G4cerr << " Can not open file " << fileName << G4endl;
00742 }
00743 #endif
00744 G4Exception( "G4ProductionCutsTable::StoreMaterialCutsCoupleInfo()",
00745 "ProcCuts102",
00746 JustWarning, "Can not open file ");
00747 return false;
00748 }
00749 G4int numberOfCouples = coupleTable.size();
00750 if (ascii) {
00752
00753 fOut << std::setw(FixedStringLengthForStore) << key << G4endl;
00754
00755
00756 fOut << numberOfCouples << G4endl;
00757 } else {
00759
00760 size_t i;
00761 for (i=0; i<FixedStringLengthForStore; ++i) temp[i] = '\0';
00762 for (i=0; i<key.length() && i<FixedStringLengthForStore-1; ++i)
00763 temp[i]=key[i];
00764 fOut.write(temp, FixedStringLengthForStore);
00765
00766
00767 fOut.write( (char*)(&numberOfCouples), sizeof (G4int));
00768 }
00769
00770
00771
00772 CoupleTableIterator cItr;
00773 for (cItr=coupleTable.begin();cItr!=coupleTable.end();cItr++){
00774 G4MaterialCutsCouple* aCouple = (*cItr);
00775 G4int index = aCouple->GetIndex();
00776
00777 G4ProductionCuts* aCut = aCouple->GetProductionCuts();
00778 G4double cutValues[NumberOfG4CutIndex];
00779 for (size_t idx=0; idx <NumberOfG4CutIndex; idx++) {
00780 cutValues[idx] = aCut->GetProductionCut(idx);
00781 }
00782
00783 G4String materialName = aCouple->GetMaterial()->GetName();
00784 G4String regionName = "NONE";
00785 if (aCouple->IsUsed()){
00786 typedef std::vector<G4Region*>::iterator regionIterator;
00787 for(regionIterator rItr=fG4RegionStore->begin();
00788 rItr!=fG4RegionStore->end();rItr++){
00789 if (IsCoupleUsedInTheRegion(aCouple, *rItr) ){
00790 regionName = (*rItr)->GetName();
00791 break;
00792 }
00793 }
00794 }
00795
00796 if (ascii) {
00798
00799 fOut << index << G4endl;
00800
00801
00802 fOut << std::setw(FixedStringLengthForStore) << materialName<< G4endl;
00803
00804
00805 fOut << std::setw(FixedStringLengthForStore) << regionName<< G4endl;
00806
00807 fOut.setf(std::ios::scientific);
00808
00809 for (size_t idx=0; idx< NumberOfG4CutIndex; idx++) {
00810 fOut << std::setw(FixedStringLengthForStore) << cutValues[idx]/(mm)
00811 << G4endl;
00812 }
00813 fOut.unsetf(std::ios::scientific);
00814
00815 } else {
00817
00818 fOut.write( (char*)(&index), sizeof (G4int));
00819
00820
00821 size_t i;
00822 for (i=0; i<FixedStringLengthForStore; ++i) temp[i] = '\0';
00823 for (i=0; i<materialName.length() && i<FixedStringLengthForStore-1; ++i) {
00824 temp[i]=materialName[i];
00825 }
00826 fOut.write(temp, FixedStringLengthForStore);
00827
00828
00829 for (i=0; i<FixedStringLengthForStore; ++i) temp[i] = '\0';
00830 for (i=0; i<regionName.length() && i<FixedStringLengthForStore-1; ++i) {
00831 temp[i]=regionName[i];
00832 }
00833 fOut.write(temp, FixedStringLengthForStore);
00834
00835
00836 for (size_t idx=0; idx< NumberOfG4CutIndex; idx++) {
00837 fOut.write( (char*)(&(cutValues[idx])), sizeof (G4double));
00838 }
00839 }
00840 }
00841
00842 fOut.close();
00843 return true;
00844 }
00845
00846
00848
00849
00850
00851 G4bool
00852 G4ProductionCutsTable::CheckMaterialCutsCoupleInfo(const G4String& directory,
00853 G4bool ascii )
00854 {
00855 const G4String fileName = directory + "/" + "couple.dat";
00856 const G4String key = "COUPLE-V3.0";
00857 std::ifstream fIn;
00858
00859
00860 if (!ascii )
00861 fIn.open(fileName,std::ios::in|std::ios::binary);
00862 else
00863 fIn.open(fileName,std::ios::in);
00864
00865
00866 if (!fIn) {
00867 #ifdef G4VERBOSE
00868 if (verboseLevel >0) {
00869 G4cerr << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo ";
00870 G4cerr << " Can not open file " << fileName << G4endl;
00871 }
00872 #endif
00873 G4Exception( "G4ProductionCutsTable::CheckMaterialCutsCoupleInfo()",
00874 "ProcCuts102",
00875 JustWarning, "Can not open file ");
00876 return false;
00877 }
00878
00879 char temp[FixedStringLengthForStore];
00880
00881
00882 G4String keyword;
00883 if (ascii) {
00884 fIn >> keyword;
00885 } else {
00886 fIn.read(temp, FixedStringLengthForStore);
00887 keyword = (const char*)(temp);
00888 }
00889 if (key!=keyword) {
00890 #ifdef G4VERBOSE
00891 if (verboseLevel >0) {
00892 G4cerr << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo ";
00893 G4cerr << " Key word in " << fileName << "= " << keyword ;
00894 G4cerr <<"( should be "<< key << ")" <<G4endl;
00895 }
00896 #endif
00897 G4Exception( "G4ProductionCutsTable::CheckMaterialCutsCoupleInfo()",
00898 "ProcCuts103",
00899 JustWarning, "Bad Data Format");
00900 fIn.close();
00901 return false;
00902 }
00903
00904
00905 G4int numberOfCouples;
00906 if (ascii) {
00907 fIn >> numberOfCouples;
00908 } else {
00909 fIn.read( (char*)(&numberOfCouples), sizeof (G4int));
00910 }
00911
00912
00913 mccConversionTable.Reset(numberOfCouples);
00914
00915
00916 for (G4int idx=0; idx<numberOfCouples; idx+=1){
00917
00918 G4int index;
00919 if (ascii) {
00920 fIn >> index;
00921 } else {
00922 fIn.read( (char*)(&index), sizeof (G4int));
00923 }
00924
00925 char mat_name[FixedStringLengthForStore];
00926 if (ascii) {
00927 fIn >> mat_name;
00928 } else {
00929 fIn.read(mat_name, FixedStringLengthForStore);
00930 }
00931
00932 char region_name[FixedStringLengthForStore];
00933 if (ascii) {
00934 fIn >> region_name;
00935 } else {
00936 fIn.read(region_name, FixedStringLengthForStore);
00937 }
00938
00939 G4double cutValues[NumberOfG4CutIndex];
00940 for (size_t i=0; i< NumberOfG4CutIndex; i++) {
00941 if (ascii) {
00942 fIn >> cutValues[i];
00943 cutValues[i] *= (mm);
00944 } else {
00945 fIn.read( (char*)(&(cutValues[i])), sizeof (G4double));
00946 }
00947 }
00948
00949
00950 CoupleTableIterator cItr;
00951 G4bool fOK = false;
00952 G4MaterialCutsCouple* aCouple =0;
00953 for (cItr=coupleTable.begin();cItr!=coupleTable.end();cItr++){
00954 aCouple = (*cItr);
00955
00956 if ( mat_name != aCouple->GetMaterial()->GetName() ) continue;
00957
00958 G4ProductionCuts* aCut = aCouple->GetProductionCuts();
00959 G4bool fRatio = true;
00960 for (size_t j=0; j< NumberOfG4CutIndex; j++) {
00961
00962 if (cutValues[j] != aCut->GetProductionCut(j)) {
00963 G4double ratio = cutValues[j]/aCut->GetProductionCut(j);
00964 fRatio = fRatio && (0.999<ratio) && (ratio<1.001) ;
00965 }
00966 }
00967 if (!fRatio) continue;
00968
00969 fOK = true;
00970 mccConversionTable.SetNewIndex(index, aCouple->GetIndex());
00971 break;
00972 }
00973
00974 #ifdef G4VERBOSE
00975
00976 if (verboseLevel >1) {
00977 if (fOK) {
00978 G4String regionname(region_name);
00979 G4Region* fRegion = 0;
00980 if ( regionname != "NONE" ) {
00981 fRegion = fG4RegionStore->GetRegion(region_name);
00982 if (fRegion==0) {
00983 G4cout << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo ";
00984 G4cout << "Region " << regionname << " is not found ";
00985 G4cout << index << ": in " << fileName << G4endl;
00986 }
00987 }
00988 if ( ( (regionname == "NONE") && (aCouple->IsUsed()) ) ||
00989 ( (fRegion !=0) && !IsCoupleUsedInTheRegion(aCouple, fRegion) ) ) {
00990 G4cout << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo ";
00991 G4cout << "A Couple is used differnt region in the current setup ";
00992 G4cout << index << ": in " << fileName << G4endl;
00993 G4cout << " material: " << mat_name ;
00994 G4cout << " region: " << region_name << G4endl;
00995 for (size_t ii=0; ii< NumberOfG4CutIndex; ii++) {
00996 G4cout << "cut[" << ii << "]=" << cutValues[ii]/mm;
00997 G4cout << " mm : ";
00998 }
00999 G4cout << G4endl;
01000 } else if ( index != aCouple->GetIndex() ) {
01001 G4cout << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo ";
01002 G4cout << "Index of couples was modified "<< G4endl;
01003 G4cout << aCouple->GetIndex() << ":" <<aCouple->GetMaterial()->GetName();
01004 G4cout <<" is defined as " ;
01005 G4cout << index << ":" << mat_name << " in " << fileName << G4endl;
01006 } else {
01007 G4cout << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo ";
01008 G4cout << index << ":" << mat_name << " in " << fileName ;
01009 G4cout << " is consistent with current setup" << G4endl;
01010 }
01011 }
01012 }
01013 if ((!fOK) && (verboseLevel >0)) {
01014 G4cout << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo ";
01015 G4cout << "Couples is not defined in the current detector setup ";
01016 G4cout << index << ": in " << fileName << G4endl;
01017 G4cout << " material: " << mat_name ;
01018 G4cout << " region: " << region_name << G4endl;
01019 for (size_t ii=0; ii< NumberOfG4CutIndex; ii++) {
01020 G4cout << "cut[" << ii << "]=" << cutValues[ii]/mm;
01021 G4cout << " mm : ";
01022 }
01023 G4cout << G4endl;
01024 }
01025 #endif
01026
01027 }
01028 fIn.close();
01029 return true;
01030 }
01031
01032
01034
01035
01036 G4bool G4ProductionCutsTable::StoreCutsInfo(const G4String& directory,
01037 G4bool ascii)
01038 {
01039 const G4String fileName = directory + "/" + "cut.dat";
01040 const G4String key = "CUT-V3.0";
01041 std::ofstream fOut;
01042 char temp[FixedStringLengthForStore];
01043
01044
01045 if (!ascii )
01046 fOut.open(fileName,std::ios::out|std::ios::binary);
01047 else
01048 fOut.open(fileName,std::ios::out);
01049
01050
01051
01052 if (!fOut) {
01053 if(verboseLevel>0) {
01054 G4cerr << "G4ProductionCutsTable::StoreCutsInfo ";
01055 G4cerr << " Can not open file " << fileName << G4endl;
01056 }
01057 G4Exception( "G4ProductionCutsTable::StoreCutsInfo()",
01058 "ProcCuts102",
01059 JustWarning, "Can not open file");
01060 return false;
01061 }
01062
01063 G4int numberOfCouples = coupleTable.size();
01064 if (ascii) {
01066
01067 fOut << key << G4endl;
01068
01069
01070 fOut << numberOfCouples << G4endl;
01071 } else {
01073
01074 size_t i;
01075 for (i=0; i<FixedStringLengthForStore; ++i) temp[i] = '\0';
01076 for (i=0; i<key.length() && i<FixedStringLengthForStore-1; ++i)
01077 temp[i]=key[i];
01078 fOut.write(temp, FixedStringLengthForStore);
01079
01080
01081 fOut.write( (char*)(&numberOfCouples), sizeof (G4int));
01082 }
01083
01084
01085 for (size_t idx=0; idx <NumberOfG4CutIndex; idx++) {
01086 const std::vector<G4double>* fRange = GetRangeCutsVector(idx);
01087 const std::vector<G4double>* fEnergy = GetEnergyCutsVector(idx);
01088 size_t i=0;
01089
01090 CoupleTableIterator cItr;
01091 for (cItr=coupleTable.begin();cItr!=coupleTable.end();cItr++, i++){
01092 if (ascii) {
01094 fOut.setf(std::ios::scientific);
01095 fOut << std::setw(20) << (*fRange)[i]/mm ;
01096 fOut << std::setw(20) << (*fEnergy)[i]/keV << G4endl;
01097 fOut.unsetf(std::ios::scientific);
01098 } else {
01100 G4double cut = (*fRange)[i];
01101 fOut.write((char*)(&cut), sizeof (G4double));
01102 cut = (*fEnergy)[i];
01103 fOut.write((char*)(&cut), sizeof (G4double));
01104 }
01105 }
01106 }
01107 fOut.close();
01108 return true;
01109 }
01110
01112
01113
01114 G4bool G4ProductionCutsTable::RetrieveCutsInfo(const G4String& directory,
01115 G4bool ascii)
01116 {
01117 const G4String fileName = directory + "/" + "cut.dat";
01118 const G4String key = "CUT-V3.0";
01119 std::ifstream fIn;
01120
01121
01122 if (!ascii )
01123 fIn.open(fileName,std::ios::in|std::ios::binary);
01124 else
01125 fIn.open(fileName,std::ios::in);
01126
01127
01128 if (!fIn) {
01129 if (verboseLevel >0) {
01130 G4cerr << "G4ProductionCutTable::RetrieveCutsInfo ";
01131 G4cerr << " Can not open file " << fileName << G4endl;
01132 }
01133 G4Exception( "G4ProductionCutsTable::RetrieveCutsInfo()",
01134 "ProcCuts102",
01135 JustWarning, "Can not open file");
01136 return false;
01137 }
01138
01139 char temp[FixedStringLengthForStore];
01140
01141
01142 G4String keyword;
01143 if (ascii) {
01144 fIn >> keyword;
01145 } else {
01146 fIn.read(temp, FixedStringLengthForStore);
01147 keyword = (const char*)(temp);
01148 }
01149 if (key!=keyword) {
01150 if (verboseLevel >0) {
01151 G4cerr << "G4ProductionCutTable::RetrieveCutsInfo ";
01152 G4cerr << " Key word in " << fileName << "= " << keyword ;
01153 G4cerr <<"( should be "<< key << ")" <<G4endl;
01154 }
01155 G4Exception( "G4ProductionCutsTable::RetrieveCutsInfo()",
01156 "ProcCuts103",
01157 JustWarning, "Bad Data Format");
01158 return false;
01159 }
01160
01161
01162 G4int numberOfCouples;
01163 if (ascii) {
01164 fIn >> numberOfCouples;
01165 if (fIn.fail()) {
01166 G4Exception( "G4ProductionCutsTable::RetrieveCutsInfo()",
01167 "ProcCuts103",
01168 JustWarning, "Bad Data Format");
01169 return false;
01170 }
01171 } else {
01172 fIn.read( (char*)(&numberOfCouples), sizeof (G4int));
01173 }
01174
01175 if (numberOfCouples > static_cast<G4int>(mccConversionTable.size()) ){
01176 G4Exception( "G4ProductionCutsTable::RetrieveCutsInfo()",
01177 "ProcCuts109", JustWarning,
01178 "Number of Couples in the file exceeds defined couples ");
01179 }
01180 numberOfCouples = mccConversionTable.size();
01181
01182 for (size_t idx=0; static_cast<G4int>(idx) <NumberOfG4CutIndex; idx++) {
01183 G4CutVectorForAParticle* fRange = rangeCutTable[idx];
01184 G4CutVectorForAParticle* fEnergy = energyCutTable[idx];
01185 fRange->clear();
01186 fEnergy->clear();
01187
01188
01189 for (size_t i=0; static_cast<G4int>(i)< numberOfCouples; i++){
01190 G4double rcut, ecut;
01191 if (ascii) {
01192 fIn >> rcut >> ecut;
01193 if (fIn.fail()) {
01194 G4Exception( "G4ProductionCutsTable::RetrieveCutsInfo()",
01195 "ProcCuts103",
01196 JustWarning, "Bad Data Format");
01197 return false;
01198 }
01199 rcut *= mm;
01200 ecut *= keV;
01201 } else {
01202 fIn.read((char*)(&rcut), sizeof (G4double));
01203 fIn.read((char*)(&ecut), sizeof (G4double));
01204 }
01205 if (!mccConversionTable.IsUsed(i)) continue;
01206 size_t new_index = mccConversionTable.GetIndex(i);
01207 (*fRange)[new_index] = rcut;
01208 (*fEnergy)[new_index] = ecut;
01209 }
01210 }
01211 return true;
01212 }
01213
01215
01216
01217 void G4ProductionCutsTable::SetVerboseLevel(G4int value)
01218 {
01219 verboseLevel = value;
01220 for (int ip=0; ip< NumberOfG4CutIndex; ip++) {
01221 if (converters[ip] !=0 ){
01222 converters[ip]->SetVerboseLevel(value);
01223 }
01224 }
01225 }
01226
01228 G4double G4ProductionCutsTable::GetMaxEnergyCut()
01229 {
01230 return G4VRangeToEnergyConverter::GetMaxEnergyCut();
01231 }
01232
01233
01235 void G4ProductionCutsTable::SetMaxEnergyCut(G4double value)
01236 {
01237 G4VRangeToEnergyConverter::SetMaxEnergyCut(value);
01238 }