Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4ProductionCutsTable.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4ProductionCutsTable.cc 78776 2014-01-23 11:26:01Z gcosmo $
28 //
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 // --------------------------------------------------------------
35 
36 #include "G4ProductionCutsTable.hh"
37 #include "G4ProductionCuts.hh"
40 #include "G4ParticleDefinition.hh"
41 #include "G4ParticleTable.hh"
42 #include "G4RegionStore.hh"
43 #include "G4LogicalVolume.hh"
44 #include "G4VPhysicalVolume.hh"
45 #include "G4RToEConvForElectron.hh"
46 #include "G4RToEConvForGamma.hh"
47 #include "G4RToEConvForPositron.hh"
48 #include "G4RToEConvForProton.hh"
49 #include "G4MaterialTable.hh"
50 #include "G4Material.hh"
51 #include "G4UnitsTable.hh"
52 
53 #include "G4Timer.hh"
54 #include "G4SystemOfUnits.hh"
55 #include "G4ios.hh"
56 #include <iomanip>
57 #include <fstream>
58 
59 
60 G4ProductionCutsTable* G4ProductionCutsTable::fG4ProductionCutsTable = 0;
61 
64 {
65  static G4ProductionCutsTable theProductionCutsTable;
66  if(!fG4ProductionCutsTable){
67  fG4ProductionCutsTable = &theProductionCutsTable;
68  }
69  return fG4ProductionCutsTable;
70 }
71 
74  : firstUse(true),verboseLevel(1),fMessenger(0)
75 {
76  for(size_t i=0;i< NumberOfG4CutIndex;i++)
77  {
78  rangeCutTable.push_back(new G4CutVectorForAParticle);
79  energyCutTable.push_back(new G4CutVectorForAParticle);
80  rangeDoubleVector[i] = 0;
81  energyDoubleVector[i] = 0;
82  converters[i] = 0;
83  }
84  fG4RegionStore = G4RegionStore::GetInstance();
85  defaultProductionCuts = new G4ProductionCuts();
86 
87  // add messenger for UI
88  fMessenger = new G4ProductionCutsTableMessenger(this);
89 }
90 
93 {
94  fMessenger=0;
95  defaultProductionCuts=0;
96  fG4RegionStore =0;
97 }
98 
101 {
102  if (defaultProductionCuts !=0) {
103  delete defaultProductionCuts;
104  defaultProductionCuts =0;
105  }
106 
107  for(CoupleTableIterator itr=coupleTable.begin();itr!=coupleTable.end();itr++){
108  delete (*itr);
109  }
110  coupleTable.clear();
111 
112  for(size_t i=0;i< NumberOfG4CutIndex;i++){
113  delete rangeCutTable[i];
114  delete energyCutTable[i];
115  delete converters[i];
116  if(rangeDoubleVector[i]!=0) delete [] rangeDoubleVector[i];
117  if(energyDoubleVector[i]!=0) delete [] energyDoubleVector[i];
118  }
119  fG4ProductionCutsTable =0;
120 
121  if (fMessenger !=0) delete fMessenger;
122  fMessenger = 0;
123 }
124 
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 
147  // Reset "used" flags of all couples
148  for(CoupleTableIterator CoupleItr=coupleTable.begin();
149  CoupleItr!=coupleTable.end();CoupleItr++)
150  {
151  (*CoupleItr)->SetUseFlag(false);
152  }
153 
154  // Update Material-Cut-Couple
155  typedef std::vector<G4Region*>::iterator regionIterator;
156  for(regionIterator rItr=fG4RegionStore->begin();
157  rItr!=fG4RegionStore->end();rItr++)
158  {
159  // Material scan is to be done only for the regions appear in the
160  // current tracking world.
161 // if((*rItr)->GetWorldPhysical()!=currentWorld) continue;
162  if((*rItr)->IsInMassGeometry() || (*rItr)->IsInParallelGeometry())
163  {
164 
165  G4ProductionCuts* fProductionCut = (*rItr)->GetProductionCuts();
166  std::vector<G4Material*>::const_iterator mItr =
167  (*rItr)->GetMaterialIterator();
168  size_t nMaterial = (*rItr)->GetNumberOfMaterials();
169  (*rItr)->ClearMap();
170 
171  for(size_t iMate=0;iMate<nMaterial;iMate++){
172  //check if this material cut couple has already been made
173  G4bool coupleAlreadyDefined = false;
174  G4MaterialCutsCouple* aCouple;
175  for(CoupleTableIterator cItr=coupleTable.begin();
176  cItr!=coupleTable.end();cItr++){
177  if( (*cItr)->GetMaterial()==(*mItr) &&
178  (*cItr)->GetProductionCuts()==fProductionCut){
179  coupleAlreadyDefined = true;
180  aCouple = *cItr;
181  break;
182  }
183  }
184 
185  // If this combination is new, cleate and register a couple
186  if(!coupleAlreadyDefined){
187  aCouple = new G4MaterialCutsCouple((*mItr),fProductionCut);
188  coupleTable.push_back(aCouple);
189  aCouple->SetIndex(coupleTable.size()-1);
190  }
191 
192  // Register this couple to the region
193  (*rItr)->RegisterMaterialCouplePair((*mItr),aCouple);
194 
195  // Set the couple to the proper logical volumes in that region
196  aCouple->SetUseFlag();
197 
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;
205 
206  ScanAndSetCouple(aLV,aCouple,aR);
207 
208  // Proceed to the next root logical volume in this region
209  rootLVItr++;
210  }
211 
212  // Proceed to next material in this region
213  mItr++;
214  }
215  }
216  }
217 
218  // Check if sizes of Range/Energy cuts tables are equal to the size of
219  // the couple table
220  // If new couples are made during the previous procedure, nCouple becomes
221  // larger then nTable
222  size_t nCouple = coupleTable.size();
223  size_t nTable = energyCutTable[0]->size();
224  G4bool newCoupleAppears = nCouple>nTable;
225  if(newCoupleAppears) {
226  for(size_t n=nCouple-nTable;n>0;n--) {
227  for(size_t nn=0;nn< NumberOfG4CutIndex;nn++){
228  rangeCutTable[nn]->push_back(-1.);
229  energyCutTable[nn]->push_back(-1.);
230  }
231  }
232  }
233 
234  // Update RangeEnergy cuts tables
235  size_t idx = 0;
236  G4Timer timer;
237  if (verboseLevel>2) {
238  timer.Start();
239  }
240  for(CoupleTableIterator cItr=coupleTable.begin();
241  cItr!=coupleTable.end();cItr++){
242  G4ProductionCuts* aCut = (*cItr)->GetProductionCuts();
243  const G4Material* aMat = (*cItr)->GetMaterial();
244  if((*cItr)->IsRecalcNeeded()) {
245  for(size_t ptcl=0;ptcl< NumberOfG4CutIndex;ptcl++){
246  G4double rCut = aCut->GetProductionCut(ptcl);
247  (*(rangeCutTable[ptcl]))[idx] = rCut;
248  // if(converters[ptcl] && (*cItr)->IsUsed())
249  if(converters[ptcl]) {
250  (*(energyCutTable[ptcl]))[idx] = converters[ptcl]->Convert(rCut,aMat);
251  } else {
252  (*(energyCutTable[ptcl]))[idx] = -1.;
253  }
254  }
255  }
256  idx++;
257  }
258  if (verboseLevel>2) {
259  timer.Stop();
260  G4cout << "G4ProductionCutsTable::UpdateCoupleTable "
261  << " elapsed time for calculation of energy cuts " << 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  }
275  }
276 
277  // Update Range/Energy cuts double vectors
278  for(size_t ix=0;ix<NumberOfG4CutIndex;ix++) {
279  for(size_t ixx=0;ixx<(*(rangeCutTable[ix])).size();ixx++) {
280  rangeDoubleVector[ix][ixx] = (*(rangeCutTable[ix]))[ixx];
281  energyDoubleVector[ix][ixx] = (*(energyCutTable[ix]))[ixx];
282  }
283  }
284 }
285 
286 
289  const G4ParticleDefinition* particle,
290  const G4Material* material,
291  G4double range
292  )
293 {
294  // This method gives energy corresponding to range value
295  // check material
296  if (material ==0) return -1.0;
297 
298  // check range
299  if (range ==0.0) return 0.0;
300  if (range <0.0) return -1.0;
301 
302  // check particle
303  G4int index = G4ProductionCuts::GetIndex(particle);
304 
305  if (index<0) {
306 #ifdef G4VERBOSE
307  if (verboseLevel >1) {
308  G4cout << "G4ProductionCutsTable::ConvertRangeToEnergy" ;
309  G4cout << particle->GetParticleName() << " has no cut value " << G4endl;
310  }
311 #endif
312  return -1.0;
313  }
314 
315  return converters[index]->Convert(range, material);
316 
317 }
318 
321 {
322  for(size_t i=0;i< NumberOfG4CutIndex;i++){
323  if (converters[i]!=0) converters[i]->Reset();
324  }
325 }
326 
327 
330 {
332 }
333 
336 {
338 }
339 
342 {
344 }
345 
346 
348 void G4ProductionCutsTable::ScanAndSetCouple(G4LogicalVolume* aLV,
349  G4MaterialCutsCouple* aCouple,
350  G4Region* aRegion)
351 {
352  //Check whether or not this logical volume belongs to the same region
353  if((aRegion!=0) && aLV->GetRegion()!=aRegion) return;
354 
355  //Check if this particular volume has a material matched to the couple
356  if(aLV->GetMaterial()==aCouple->GetMaterial()) {
357  aLV->SetMaterialCutsCouple(aCouple);
358  }
359 
360  size_t noDaughters = aLV->GetNoDaughters();
361  if(noDaughters==0) return;
362 
363  //Loop over daughters with same region
364  for(size_t i=0;i<noDaughters;i++){
365  G4LogicalVolume* daughterLVol = aLV->GetDaughter(i)->GetLogicalVolume();
366  ScanAndSetCouple(daughterLVol,aCouple,aRegion);
367  }
368 }
369 
372 {
373  G4cout << G4endl;
374  G4cout << "========= Table of registered couples =============================="
375  << G4endl;
376  for(CoupleTableIterator cItr=coupleTable.begin();
377  cItr!=coupleTable.end();cItr++) {
378  G4MaterialCutsCouple* aCouple = (*cItr);
379  G4ProductionCuts* aCut = aCouple->GetProductionCuts();
380  G4cout << G4endl;
381  G4cout << "Index : " << aCouple->GetIndex()
382  << " used in the geometry : ";
383  if(aCouple->IsUsed()) G4cout << "Yes";
384  else G4cout << "No ";
388  G4cout << G4endl;
389  G4cout << " Material : " << aCouple->GetMaterial()->GetName() << G4endl;
390  G4cout << " Range cuts : "
391  << " gamma " << G4BestUnit(aCut->GetProductionCut("gamma"),"Length")
392  << " e- " << G4BestUnit(aCut->GetProductionCut("e-"),"Length")
393  << " e+ " << G4BestUnit(aCut->GetProductionCut("e+"),"Length")
394  << " proton " << G4BestUnit(aCut->GetProductionCut("proton"),"Length")
395  << G4endl;
396  G4cout << " Energy thresholds : " ;
400  G4cout << " gamma " << G4BestUnit((*(energyCutTable[0]))[aCouple->GetIndex()],"Energy")
401  << " e- " << G4BestUnit((*(energyCutTable[1]))[aCouple->GetIndex()],"Energy")
402  << " e+ " << G4BestUnit((*(energyCutTable[2]))[aCouple->GetIndex()],"Energy")
403  << " proton " << G4BestUnit((*(energyCutTable[3]))[aCouple->GetIndex()],"Energy");
405  G4cout << G4endl;
406 
407  if(aCouple->IsUsed()) {
408  G4cout << " Region(s) which use this couple : " << G4endl;
409  typedef std::vector<G4Region*>::iterator regionIterator;
410  for(regionIterator rItr=fG4RegionStore->begin();
411  rItr!=fG4RegionStore->end();rItr++) {
412  if (IsCoupleUsedInTheRegion(aCouple, *rItr) ){
413  G4cout << " " << (*rItr)->GetName() << G4endl;
414  }
415  }
416  }
417  }
418  G4cout << G4endl;
419  G4cout << "====================================================================" << G4endl;
420  G4cout << G4endl;
421 }
422 
423 
425 // Store cuts and material information in files under the specified directory.
427  G4bool ascii)
428 {
429  if (!StoreMaterialInfo(dir, ascii)) return false;
430  if (!StoreMaterialCutsCoupleInfo(dir, ascii)) return false;
431  if (!StoreCutsInfo(dir, ascii)) return false;
432 
433 #ifdef G4VERBOSE
434  if (verboseLevel >2) {
435  G4cout << "G4ProductionCutsTable::StoreCutsTable " ;
436  G4cout << " Material/Cuts information have been succesfully stored ";
437  if (ascii) {
438  G4cout << " in Ascii mode ";
439  }else{
440  G4cout << " in Binary mode ";
441  }
442  G4cout << " under " << dir << G4endl;
443  }
444 #endif
445  return true;
446 }
447 
450  G4bool ascii)
451 {
452  if (!CheckForRetrieveCutsTable(dir, ascii)) return false;
453  if (!RetrieveCutsInfo(dir, ascii)) return false;
454 #ifdef G4VERBOSE
455  if (verboseLevel >2) {
456  G4cout << "G4ProductionCutsTable::RetrieveCutsTable " ;
457  G4cout << " Material/Cuts information have been succesfully retreived ";
458  if (ascii) {
459  G4cout << " in Ascii mode ";
460  }else{
461  G4cout << " in Binary mode ";
462  }
463  G4cout << " under " << dir << G4endl;
464  }
465 #endif
466  return true;
467 }
468 
470 // check stored material and cut values are consistent
471 // with the current detector setup.
472 //
473 G4bool
475  G4bool ascii)
476 {
477  G4cerr << "G4ProductionCutsTable::CheckForRetrieveCutsTable!!"<< G4endl;
478  // isNeedForRestoreCoupleInfo = false;
479  if (!CheckMaterialInfo(directory, ascii)) return false;
480  if (verboseLevel >2) {
481  G4cerr << "G4ProductionCutsTable::CheckMaterialInfo passed !!"<< G4endl;
482  }
483  if (!CheckMaterialCutsCoupleInfo(directory, ascii)) return false;
484  if (verboseLevel >2) {
485  G4cerr << "G4ProductionCutsTable::CheckMaterialCutsCoupleInfo passed !!"<< G4endl;
486  }
487  return true;
488 }
489 
491 // Store material information in files under the specified directory.
492 //
494  G4bool ascii)
495 {
496  const G4String fileName = directory + "/" + "material.dat";
497  const G4String key = "MATERIAL-V3.0";
498  std::ofstream fOut;
499 
500  // open output file //
501  if (!ascii )
502  fOut.open(fileName,std::ios::out|std::ios::binary);
503  else
504  fOut.open(fileName,std::ios::out);
505 
506 
507  // check if the file has been opened successfully
508  if (!fOut) {
509 #ifdef G4VERBOSE
510  if (verboseLevel>0) {
511  G4cerr << "G4ProductionCutsTable::StoreMaterialInfo ";
512  G4cerr << " Can not open file " << fileName << G4endl;
513  }
514 #endif
515  G4Exception( "G4ProductionCutsTable::StoreMaterialInfo()",
516  "ProcCuts102",
517  JustWarning, "Can not open file ");
518  return false;
519  }
520 
521  const G4MaterialTable* matTable = G4Material::GetMaterialTable();
522  // number of materials in the table
523  G4int numberOfMaterial = matTable->size();
524 
525  if (ascii) {
527  // key word
528  fOut << key << G4endl;
529 
530  // number of materials in the table
531  fOut << numberOfMaterial << G4endl;
532 
533  fOut.setf(std::ios::scientific);
534 
535  // material name and density
536  for (size_t idx=0; static_cast<G4int>(idx)<numberOfMaterial; ++idx){
537  fOut << std::setw(FixedStringLengthForStore)
538  << ((*matTable)[idx])->GetName();
539  fOut << std::setw(FixedStringLengthForStore)
540  << ((*matTable)[idx])->GetDensity()/(g/cm3) << G4endl;
541  }
542 
543  fOut.unsetf(std::ios::scientific);
544 
545  } else {
547  char temp[FixedStringLengthForStore];
548  size_t i;
549 
550  // key word
551  for (i=0; i<FixedStringLengthForStore; ++i){
552  temp[i] = '\0';
553  }
554  for (i=0; i<key.length() && i<FixedStringLengthForStore-1; ++i){
555  temp[i]=key[i];
556  }
557  fOut.write(temp, FixedStringLengthForStore);
558 
559  // number of materials in the table
560  fOut.write( (char*)(&numberOfMaterial), sizeof (G4int));
561 
562  // material name and density
563  for (size_t imat=0; static_cast<G4int>(imat)<numberOfMaterial; ++imat){
564  G4String name = ((*matTable)[imat])->GetName();
565  G4double density = ((*matTable)[imat])->GetDensity();
566  for (i=0; i<FixedStringLengthForStore; ++i) temp[i] = '\0';
567  for (i=0; i<name.length() && i<FixedStringLengthForStore-1; ++i)
568  temp[i]=name[i];
569  fOut.write(temp, FixedStringLengthForStore);
570  fOut.write( (char*)(&density), sizeof (G4double));
571  }
572  }
573 
574  fOut.close();
575  return true;
576 }
577 
579 // check stored material is consistent with the current detector setup.
581  G4bool ascii)
582 {
583  const G4String fileName = directory + "/" + "material.dat";
584  const G4String key = "MATERIAL-V3.0";
585  std::ifstream fIn;
586 
587  // open input file //
588  if (!ascii )
589  fIn.open(fileName,std::ios::in|std::ios::binary);
590  else
591  fIn.open(fileName,std::ios::in);
592 
593  // check if the file has been opened successfully
594  if (!fIn) {
595 #ifdef G4VERBOSE
596  if (verboseLevel >0) {
597  G4cerr << "G4ProductionCutsTable::CheckMaterialInfo ";
598  G4cerr << " Can not open file " << fileName << G4endl;
599  }
600 #endif
601  G4Exception( "G4ProductionCutsTable::CheckMaterialInfo()",
602  "ProcCuts102",
603  JustWarning, "Can not open file ");
604  return false;
605  }
606 
607  char temp[FixedStringLengthForStore];
608 
609  // key word
610  G4String keyword;
611  if (ascii) {
612  fIn >> keyword;
613  } else {
614  fIn.read(temp, FixedStringLengthForStore);
615  keyword = (const char*)(temp);
616  }
617  if (key!=keyword) {
618 #ifdef G4VERBOSE
619  if (verboseLevel >0) {
620  G4cerr << "G4ProductionCutsTable::CheckMaterialInfo ";
621  G4cerr << " Key word in " << fileName << "= " << keyword ;
622  G4cerr <<"( should be "<< key << ")" <<G4endl;
623  }
624 #endif
625  G4Exception( "G4ProductionCutsTable::CheckMaterialInfo()",
626  "ProcCuts103",
627  JustWarning, "Bad Data Format");
628  return false;
629  }
630 
631  // number of materials in the table
632  G4int nmat;
633  if (ascii) {
634  fIn >> nmat;
635  } else {
636  fIn.read( (char*)(&nmat), sizeof (G4int));
637  }
638  if ((nmat<=0) || (nmat >100000)){
639  G4Exception( "G4ProductionCutsTable::CheckMaterialInfo()",
640  "ProcCuts108",JustWarning,
641  "Number of materials is less than zero or too big");
642  return false;
643  }
644 
645  // list of material
646  for (G4int idx=0; idx<nmat ; ++idx){
647  // check eof
648  if(fIn.eof()) {
649 #ifdef G4VERBOSE
650  if (verboseLevel >0) {
651  G4cout << "G4ProductionCutsTable::CheckMaterialInfo ";
652  G4cout << " encountered End of File " ;
653  G4cout << " at " << idx+1 << "th material "<< G4endl;
654  }
655 #endif
656  fIn.close();
657  return false;
658  }
659 
660  // check material name and density
661  char name[FixedStringLengthForStore];
662  double density;
663  if (ascii) {
664  fIn >> name >> density;
665  density *= (g/cm3);
666 
667  } else {
668  fIn.read(name, FixedStringLengthForStore);
669  fIn.read((char*)(&density), sizeof (G4double));
670  }
671  if (fIn.fail()) {
672 #ifdef G4VERBOSE
673  if (verboseLevel >0) {
674  G4cerr << "G4ProductionCutsTable::CheckMaterialInfo ";
675  G4cerr << " Bad data format ";
676  G4cerr << " at " << idx+1 << "th material "<< G4endl;
677  }
678 #endif
679  G4Exception( "G4ProductionCutsTable::CheckMaterialInfo()",
680  "ProcCuts103",
681  JustWarning, "Bad Data Format");
682  fIn.close();
683  return false;
684  }
685 
686  G4Material* aMaterial = G4Material::GetMaterial(name);
687  if (aMaterial ==0 ) continue;
688 
689  G4double ratio = std::fabs(density/aMaterial->GetDensity() );
690  if ((0.999>ratio) || (ratio>1.001) ){
691 #ifdef G4VERBOSE
692  if (verboseLevel >0) {
693  G4cerr << "G4ProductionCutsTable::CheckMaterialInfo ";
694  G4cerr << " Inconsistent material density" << G4endl;;
695  G4cerr << " at " << idx+1 << "th material "<< G4endl;
696  G4cerr << "Name: " << name << G4endl;
697  G4cerr << "Density:" << std::setiosflags(std::ios::scientific) << density / (g/cm3) ;
698  G4cerr << "(should be " << aMaterial->GetDensity()/(g/cm3)<< ")" << " [g/cm3]"<< G4endl;
699  G4cerr << std::resetiosflags(std::ios::scientific);
700  }
701 #endif
702  G4Exception( "G4ProductionCutsTable::CheckMaterialInfo()",
703  "ProcCuts104",
704  JustWarning, "Inconsitent matrial density");
705  fIn.close();
706  return false;
707  }
708 
709  }
710 
711  fIn.close();
712  return true;
713 
714 }
715 
716 
718 // Store materialCutsCouple information in files under the specified directory.
719 //
720 G4bool
722  G4bool ascii)
723 {
724  const G4String fileName = directory + "/" + "couple.dat";
725  const G4String key = "COUPLE-V3.0";
726  std::ofstream fOut;
727  char temp[FixedStringLengthForStore];
728 
729  // open output file //
730  if (!ascii )
731  fOut.open(fileName,std::ios::out|std::ios::binary);
732  else
733  fOut.open(fileName,std::ios::out);
734 
735 
736  // check if the file has been opened successfully
737  if (!fOut) {
738 #ifdef G4VERBOSE
739  if (verboseLevel >0) {
740  G4cerr << "G4ProductionCutsTable::StoreMaterialCutsCoupleInfo ";
741  G4cerr << " Can not open file " << fileName << G4endl;
742  }
743 #endif
744  G4Exception( "G4ProductionCutsTable::StoreMaterialCutsCoupleInfo()",
745  "ProcCuts102",
746  JustWarning, "Can not open file ");
747  return false;
748  }
749  G4int numberOfCouples = coupleTable.size();
750  if (ascii) {
752  // key word
753  fOut << std::setw(FixedStringLengthForStore) << key << G4endl;
754 
755  // number of couples in the table
756  fOut << numberOfCouples << G4endl;
757  } else {
759  // key word
760  size_t i;
761  for (i=0; i<FixedStringLengthForStore; ++i) temp[i] = '\0';
762  for (i=0; i<key.length() && i<FixedStringLengthForStore-1; ++i)
763  temp[i]=key[i];
764  fOut.write(temp, FixedStringLengthForStore);
765 
766  // number of couples in the table
767  fOut.write( (char*)(&numberOfCouples), sizeof (G4int));
768  }
769 
770 
771  // Loop over all couples
772  CoupleTableIterator cItr;
773  for (cItr=coupleTable.begin();cItr!=coupleTable.end();cItr++){
774  G4MaterialCutsCouple* aCouple = (*cItr);
775  G4int index = aCouple->GetIndex();
776  // cut value
777  G4ProductionCuts* aCut = aCouple->GetProductionCuts();
778  G4double cutValues[NumberOfG4CutIndex];
779  for (size_t idx=0; idx <NumberOfG4CutIndex; idx++) {
780  cutValues[idx] = aCut->GetProductionCut(idx);
781  }
782  // material/region info
783  G4String materialName = aCouple->GetMaterial()->GetName();
784  G4String regionName = "NONE";
785  if (aCouple->IsUsed()){
786  typedef std::vector<G4Region*>::iterator regionIterator;
787  for(regionIterator rItr=fG4RegionStore->begin();
788  rItr!=fG4RegionStore->end();rItr++){
789  if (IsCoupleUsedInTheRegion(aCouple, *rItr) ){
790  regionName = (*rItr)->GetName();
791  break;
792  }
793  }
794  }
795 
796  if (ascii) {
798  // index number
799  fOut << index << G4endl;
800 
801  // material name
802  fOut << std::setw(FixedStringLengthForStore) << materialName<< G4endl;
803 
804  // region name
805  fOut << std::setw(FixedStringLengthForStore) << regionName<< G4endl;
806 
807  fOut.setf(std::ios::scientific);
808  // cut values
809  for (size_t idx=0; idx< NumberOfG4CutIndex; idx++) {
810  fOut << std::setw(FixedStringLengthForStore) << cutValues[idx]/(mm)
811  << G4endl;
812  }
813  fOut.unsetf(std::ios::scientific);
814 
815  } else {
817  // index
818  fOut.write( (char*)(&index), sizeof (G4int));
819 
820  // material name
821  size_t i;
822  for (i=0; i<FixedStringLengthForStore; ++i) temp[i] = '\0';
823  for (i=0; i<materialName.length() && i<FixedStringLengthForStore-1; ++i) {
824  temp[i]=materialName[i];
825  }
826  fOut.write(temp, FixedStringLengthForStore);
827 
828  // region name
829  for (i=0; i<FixedStringLengthForStore; ++i) temp[i] = '\0';
830  for (i=0; i<regionName.length() && i<FixedStringLengthForStore-1; ++i) {
831  temp[i]=regionName[i];
832  }
833  fOut.write(temp, FixedStringLengthForStore);
834 
835  // cut values
836  for (size_t idx=0; idx< NumberOfG4CutIndex; idx++) {
837  fOut.write( (char*)(&(cutValues[idx])), sizeof (G4double));
838  }
839  }
840  }
841 
842  fOut.close();
843  return true;
844 }
845 
846 
848 // check stored materialCutsCouple is consistent
849 // with the current detector setup.
850 //
851 G4bool
853  G4bool ascii )
854 {
855  const G4String fileName = directory + "/" + "couple.dat";
856  const G4String key = "COUPLE-V3.0";
857  std::ifstream fIn;
858 
859  // open input file //
860  if (!ascii )
861  fIn.open(fileName,std::ios::in|std::ios::binary);
862  else
863  fIn.open(fileName,std::ios::in);
864 
865  // check if the file has been opened successfully
866  if (!fIn) {
867 #ifdef G4VERBOSE
868  if (verboseLevel >0) {
869  G4cerr << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo ";
870  G4cerr << " Can not open file " << fileName << G4endl;
871  }
872 #endif
873  G4Exception( "G4ProductionCutsTable::CheckMaterialCutsCoupleInfo()",
874  "ProcCuts102",
875  JustWarning, "Can not open file ");
876  return false;
877  }
878 
879  char temp[FixedStringLengthForStore];
880 
881  // key word
882  G4String keyword;
883  if (ascii) {
884  fIn >> keyword;
885  } else {
886  fIn.read(temp, FixedStringLengthForStore);
887  keyword = (const char*)(temp);
888  }
889  if (key!=keyword) {
890 #ifdef G4VERBOSE
891  if (verboseLevel >0) {
892  G4cerr << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo ";
893  G4cerr << " Key word in " << fileName << "= " << keyword ;
894  G4cerr <<"( should be "<< key << ")" <<G4endl;
895  }
896 #endif
897  G4Exception( "G4ProductionCutsTable::CheckMaterialCutsCoupleInfo()",
898  "ProcCuts103",
899  JustWarning, "Bad Data Format");
900  fIn.close();
901  return false;
902  }
903 
904  // numberOfCouples
905  G4int numberOfCouples;
906  if (ascii) {
907  fIn >> numberOfCouples;
908  } else {
909  fIn.read( (char*)(&numberOfCouples), sizeof (G4int));
910  }
911 
912  // Reset MCCIndexConversionTable
913  mccConversionTable.Reset(numberOfCouples);
914 
915  // Read in couple information
916  for (G4int idx=0; idx<numberOfCouples; idx+=1){
917  // read in index
918  G4int index;
919  if (ascii) {
920  fIn >> index;
921  } else {
922  fIn.read( (char*)(&index), sizeof (G4int));
923  }
924  // read in index material name
925  char mat_name[FixedStringLengthForStore];
926  if (ascii) {
927  fIn >> mat_name;
928  } else {
929  fIn.read(mat_name, FixedStringLengthForStore);
930  }
931  // read in index and region name
932  char region_name[FixedStringLengthForStore];
933  if (ascii) {
934  fIn >> region_name;
935  } else {
936  fIn.read(region_name, FixedStringLengthForStore);
937  }
938  // cut value
939  G4double cutValues[NumberOfG4CutIndex];
940  for (size_t i=0; i< NumberOfG4CutIndex; i++) {
941  if (ascii) {
942  fIn >> cutValues[i];
943  cutValues[i] *= (mm);
944  } else {
945  fIn.read( (char*)(&(cutValues[i])), sizeof (G4double));
946  }
947  }
948 
949  // Loop over all couples
950  CoupleTableIterator cItr;
951  G4bool fOK = false;
952  G4MaterialCutsCouple* aCouple =0;
953  for (cItr=coupleTable.begin();cItr!=coupleTable.end();cItr++){
954  aCouple = (*cItr);
955  // check material name
956  if ( mat_name != aCouple->GetMaterial()->GetName() ) continue;
957  // check cut values
958  G4ProductionCuts* aCut = aCouple->GetProductionCuts();
959  G4bool fRatio = true;
960  for (size_t j=0; j< NumberOfG4CutIndex; j++) {
961  // check ratio only if values are not the same
962  if (cutValues[j] != aCut->GetProductionCut(j)) {
963  G4double ratio = cutValues[j]/aCut->GetProductionCut(j);
964  fRatio = fRatio && (0.999<ratio) && (ratio<1.001) ;
965  }
966  }
967  if (!fRatio) continue;
968  // MCC matched
969  fOK = true;
970  mccConversionTable.SetNewIndex(index, aCouple->GetIndex());
971  break;
972  }
973 
974  if (fOK) {
975 #ifdef G4VERBOSE
976  // debug information
977  if (verboseLevel >1) {
978  G4String regionname(region_name);
979  G4Region* fRegion = 0;
980  if ( regionname != "NONE" ) {
981  fRegion = fG4RegionStore->GetRegion(region_name);
982  if (fRegion==0) {
983  G4cout << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo ";
984  G4cout << "Region " << regionname << " is not found ";
985  G4cout << index << ": in " << fileName << G4endl;
986  }
987  }
988  if ( ( (regionname == "NONE") && (aCouple->IsUsed()) ) ||
989  ( (fRegion !=0) && !IsCoupleUsedInTheRegion(aCouple, fRegion) ) ) {
990  G4cout << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo ";
991  G4cout << "A Couple is used differnt region in the current setup ";
992  G4cout << index << ": in " << fileName << G4endl;
993  G4cout << " material: " << mat_name ;
994  G4cout << " region: " << region_name << G4endl;
995  for (size_t ii=0; ii< NumberOfG4CutIndex; ii++) {
996  G4cout << "cut[" << ii << "]=" << cutValues[ii]/mm;
997  G4cout << " mm : ";
998  }
999  G4cout << G4endl;
1000  } else if ( index != aCouple->GetIndex() ) {
1001  G4cout << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo ";
1002  G4cout << "Index of couples was modified "<< G4endl;
1003  G4cout << aCouple->GetIndex() << ":" <<aCouple->GetMaterial()->GetName();
1004  G4cout <<" is defined as " ;
1005  G4cout << index << ":" << mat_name << " in " << fileName << G4endl;
1006  } else {
1007  G4cout << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo ";
1008  G4cout << index << ":" << mat_name << " in " << fileName ;
1009  G4cout << " is consistent with current setup" << G4endl;
1010  }
1011  }
1012 #endif
1013  } else {
1014 #ifdef G4VERBOSE
1015  if (verboseLevel >0) {
1016  G4cout << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo ";
1017  G4cout << "Couples is not defined in the current detector setup ";
1018  G4cout << index << ": in " << fileName << G4endl;
1019  G4cout << " material: " << mat_name ;
1020  G4cout << " region: " << region_name << G4endl;
1021  for (size_t ii=0; ii< NumberOfG4CutIndex; ii++) {
1022  G4cout << "cut[" << ii << "]=" << cutValues[ii]/mm;
1023  G4cout << " mm : ";
1024  }
1025  G4cout << G4endl;
1026  }
1027 #endif
1028  }
1029 
1030  }
1031  fIn.close();
1032  return true;
1033 }
1034 
1035 
1037 // Store cut values information in files under the specified directory.
1038 //
1040  G4bool ascii)
1041 {
1042  const G4String fileName = directory + "/" + "cut.dat";
1043  const G4String key = "CUT-V3.0";
1044  std::ofstream fOut;
1045  char temp[FixedStringLengthForStore];
1046 
1047  // open output file //
1048  if (!ascii )
1049  fOut.open(fileName,std::ios::out|std::ios::binary);
1050  else
1051  fOut.open(fileName,std::ios::out);
1052 
1053 
1054  // check if the file has been opened successfully
1055  if (!fOut) {
1056  if(verboseLevel>0) {
1057  G4cerr << "G4ProductionCutsTable::StoreCutsInfo ";
1058  G4cerr << " Can not open file " << fileName << G4endl;
1059  }
1060  G4Exception( "G4ProductionCutsTable::StoreCutsInfo()",
1061  "ProcCuts102",
1062  JustWarning, "Can not open file");
1063  return false;
1064  }
1065 
1066  G4int numberOfCouples = coupleTable.size();
1067  if (ascii) {
1069  // key word
1070  fOut << key << G4endl;
1071 
1072  // number of couples in the table
1073  fOut << numberOfCouples << G4endl;
1074  } else {
1076  // key word
1077  size_t i;
1078  for (i=0; i<FixedStringLengthForStore; ++i) temp[i] = '\0';
1079  for (i=0; i<key.length() && i<FixedStringLengthForStore-1; ++i)
1080  temp[i]=key[i];
1081  fOut.write(temp, FixedStringLengthForStore);
1082 
1083  // number of couples in the table
1084  fOut.write( (char*)(&numberOfCouples), sizeof (G4int));
1085  }
1086 
1087 
1088  for (size_t idx=0; idx <NumberOfG4CutIndex; idx++) {
1089  const std::vector<G4double>* fRange = GetRangeCutsVector(idx);
1090  const std::vector<G4double>* fEnergy = GetEnergyCutsVector(idx);
1091  size_t i=0;
1092  // Loop over all couples
1093  CoupleTableIterator cItr;
1094  for (cItr=coupleTable.begin();cItr!=coupleTable.end();cItr++, i++){
1095  if (ascii) {
1097  fOut.setf(std::ios::scientific);
1098  fOut << std::setw(20) << (*fRange)[i]/mm ;
1099  fOut << std::setw(20) << (*fEnergy)[i]/keV << G4endl;
1100  fOut.unsetf(std::ios::scientific);
1101  } else {
1103  G4double cut = (*fRange)[i];
1104  fOut.write((char*)(&cut), sizeof (G4double));
1105  cut = (*fEnergy)[i];
1106  fOut.write((char*)(&cut), sizeof (G4double));
1107  }
1108  }
1109  }
1110  fOut.close();
1111  return true;
1112 }
1113 
1115 // Retrieve cut values information in files under the specified directory.
1116 //
1118  G4bool ascii)
1119 {
1120  const G4String fileName = directory + "/" + "cut.dat";
1121  const G4String key = "CUT-V3.0";
1122  std::ifstream fIn;
1123 
1124  // open input file //
1125  if (!ascii )
1126  fIn.open(fileName,std::ios::in|std::ios::binary);
1127  else
1128  fIn.open(fileName,std::ios::in);
1129 
1130  // check if the file has been opened successfully
1131  if (!fIn) {
1132  if (verboseLevel >0) {
1133  G4cerr << "G4ProductionCutTable::RetrieveCutsInfo ";
1134  G4cerr << " Can not open file " << fileName << G4endl;
1135  }
1136  G4Exception( "G4ProductionCutsTable::RetrieveCutsInfo()",
1137  "ProcCuts102",
1138  JustWarning, "Can not open file");
1139  return false;
1140  }
1141 
1142  char temp[FixedStringLengthForStore];
1143 
1144  // key word
1145  G4String keyword;
1146  if (ascii) {
1147  fIn >> keyword;
1148  } else {
1149  fIn.read(temp, FixedStringLengthForStore);
1150  keyword = (const char*)(temp);
1151  }
1152  if (key!=keyword) {
1153  if (verboseLevel >0) {
1154  G4cerr << "G4ProductionCutTable::RetrieveCutsInfo ";
1155  G4cerr << " Key word in " << fileName << "= " << keyword ;
1156  G4cerr <<"( should be "<< key << ")" <<G4endl;
1157  }
1158  G4Exception( "G4ProductionCutsTable::RetrieveCutsInfo()",
1159  "ProcCuts103",
1160  JustWarning, "Bad Data Format");
1161  return false;
1162  }
1163 
1164  // numberOfCouples
1165  G4int numberOfCouples;
1166  if (ascii) {
1167  fIn >> numberOfCouples;
1168  if (fIn.fail()) {
1169  G4Exception( "G4ProductionCutsTable::RetrieveCutsInfo()",
1170  "ProcCuts103",
1171  JustWarning, "Bad Data Format");
1172  return false;
1173  }
1174  } else {
1175  fIn.read( (char*)(&numberOfCouples), sizeof (G4int));
1176  }
1177 
1178  if (numberOfCouples > static_cast<G4int>(mccConversionTable.size()) ){
1179  G4Exception( "G4ProductionCutsTable::RetrieveCutsInfo()",
1180  "ProcCuts109", JustWarning,
1181  "Number of Couples in the file exceeds defined couples ");
1182  }
1183  numberOfCouples = mccConversionTable.size();
1184 
1185  for (size_t idx=0; static_cast<G4int>(idx) <NumberOfG4CutIndex; idx++) {
1186  G4CutVectorForAParticle* fRange = rangeCutTable[idx];
1187  G4CutVectorForAParticle* fEnergy = energyCutTable[idx];
1188  fRange->clear();
1189  fEnergy->clear();
1190 
1191  // Loop over all couples
1192  for (size_t i=0; static_cast<G4int>(i)< numberOfCouples; i++){
1193  G4double rcut, ecut;
1194  if (ascii) {
1195  fIn >> rcut >> ecut;
1196  if (fIn.fail()) {
1197  G4Exception( "G4ProductionCutsTable::RetrieveCutsInfo()",
1198  "ProcCuts103",
1199  JustWarning, "Bad Data Format");
1200  return false;
1201  }
1202  rcut *= mm;
1203  ecut *= keV;
1204  } else {
1205  fIn.read((char*)(&rcut), sizeof (G4double));
1206  fIn.read((char*)(&ecut), sizeof (G4double));
1207  }
1208  if (!mccConversionTable.IsUsed(i)) continue;
1209  size_t new_index = mccConversionTable.GetIndex(i);
1210  (*fRange)[new_index] = rcut;
1211  (*fEnergy)[new_index] = ecut;
1212  }
1213  }
1214  return true;
1215 }
1216 
1218 // Set Verbose Level
1219 // set same verbosity to all registered RangeToEnergyConverters
1221 {
1222  verboseLevel = value;
1223  for (int ip=0; ip< NumberOfG4CutIndex; ip++) {
1224  if (converters[ip] !=0 ){
1225  converters[ip]->SetVerboseLevel(value);
1226  }
1227  }
1228 }
1229 
1232 {
1234 }
1235 
1236 
1239 {
1241 }
const XML_Char * name
Definition: expat.h:151
virtual G4bool StoreMaterialInfo(const G4String &directory, G4bool ascii=false)
G4Region * GetRegion(const G4String &name, G4bool verbose=true) const
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
G4double GetHighEdgeEnergy() const
static G4int GetIndex(const G4String &name)
static constexpr double mm
Definition: G4SIunits.hh:115
G4Material * GetMaterial() const
G4int GetIndex(size_t index) const
void SetEnergyRange(G4double lowedge, G4double highedge)
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:602
G4bool CheckForRetrieveCutsTable(const G4String &directory, G4bool ascii=false)
G4double GetProductionCut(G4int index) const
const G4String & GetName() const
Definition: G4Material.hh:178
void SetNewIndex(size_t index, size_t new_value)
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:587
std::vector< G4Material * > G4MaterialTable
G4double GetDensity() const
Definition: G4Material.hh:180
G4VPhysicalVolume * GetDaughter(const G4int i) const
virtual G4bool CheckMaterialCutsCoupleInfo(const G4String &directory, G4bool ascii=false)
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
G4Region * GetRegion() const
static constexpr double g
Definition: G4SIunits.hh:183
int G4int
Definition: G4Types.hh:78
G4bool IsUsed(size_t index) const
const G4String & GetParticleName() const
static G4RegionStore * GetInstance()
G4GLOB_DLL std::ostream G4cout
virtual G4double Convert(G4double rangeCut, const G4Material *material)
const XML_Char int const XML_Char * value
Definition: expat.h:331
virtual G4bool RetrieveCutsInfo(const G4String &directory, G4bool ascii=false)
void SetVerboseLevel(G4int value)
G4bool RetrieveCutsTable(const G4String &directory, G4bool ascii=false)
bool G4bool
Definition: G4Types.hh:79
static void SetMaxEnergyCut(G4double value)
void SetMaterialCutsCouple(G4MaterialCutsCouple *cuts)
static void SetEnergyRange(G4double lowedge, G4double highedge)
const std::vector< G4double > * GetRangeCutsVector(size_t pcIdx) const
static constexpr double cm3
Definition: G4SIunits.hh:121
G4int GetNoDaughters() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
const std::vector< G4double > & GetProductionCuts() const
G4double GetLowEdgeEnergy() const
static G4ProductionCutsTable * GetProductionCutsTable()
G4LogicalVolume * GetLogicalVolume() const
static G4ParticleTable * GetParticleTable()
void Stop()
virtual G4bool StoreMaterialCutsCoupleInfo(const G4String &directory, G4bool ascii=false)
void SetUseFlag(G4bool flg=true)
G4bool StoreCutsTable(const G4String &directory, G4bool ascii=false)
#define G4endl
Definition: G4ios.hh:61
void UpdateCoupleTable(G4VPhysicalVolume *currentWorld)
G4double ConvertRangeToEnergy(const G4ParticleDefinition *particle, const G4Material *material, G4double range)
void Start()
double G4double
Definition: G4Types.hh:76
virtual G4bool CheckMaterialInfo(const G4String &directory, G4bool ascii=false)
G4ProductionCuts * GetProductionCuts() const
static constexpr double keV
Definition: G4SIunits.hh:216
void SetMaxEnergyCut(G4double value)
virtual G4bool StoreCutsInfo(const G4String &directory, G4bool ascii=false)
const G4Material * GetMaterial() const
G4GLOB_DLL std::ostream G4cerr