source: trunk/source/processes/cuts/src/G4ProductionCutsTable.cc@ 1272

Last change on this file since 1272 was 1228, checked in by garnier, 16 years ago

update geant4.9.3 tag
File size: 36.7 KB
RevLine 
[819]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//
[1196]27// $Id: G4ProductionCutsTable.cc,v 1.25 2009/11/11 03:20:22 kurasige Exp $
[1228]28// GEANT4 tag $Name: geant4-09-03 $
[819]29//
30//
31// --------------------------------------------------------------
32// GEANT 4 class implementation file/ History:
33// 06/Oct. 2002, M.Asai : First implementation
[961]34// 02/Mar. 2008, H.Kurashige : Add messenger
[819]35// --------------------------------------------------------------
36
37#include "G4ProductionCutsTable.hh"
38#include "G4ProductionCuts.hh"
39#include "G4MCCIndexConversionTable.hh"
[961]40#include "G4ProductionCutsTableMessenger.hh"
[819]41#include "G4ParticleDefinition.hh"
42#include "G4ParticleTable.hh"
43#include "G4RegionStore.hh"
44#include "G4LogicalVolume.hh"
45#include "G4VPhysicalVolume.hh"
46#include "G4RToEConvForElectron.hh"
47#include "G4RToEConvForGamma.hh"
48#include "G4RToEConvForPositron.hh"
[1196]49#include "G4RToEConvForProton.hh"
[819]50#include "G4MaterialTable.hh"
51#include "G4Material.hh"
52#include "G4UnitsTable.hh"
53
[1196]54#include "G4Timer.hh"
55
[819]56#include "G4ios.hh"
57#include <iomanip>
58#include <fstream>
59
60
61G4ProductionCutsTable* G4ProductionCutsTable::fG4ProductionCutsTable = 0;
62
[1196]63/////////////////////////////////////////////////////////////
[819]64G4ProductionCutsTable* G4ProductionCutsTable::GetProductionCutsTable()
65{
66 static G4ProductionCutsTable theProductionCutsTable;
67 if(!fG4ProductionCutsTable){
68 fG4ProductionCutsTable = &theProductionCutsTable;
69 }
70 return fG4ProductionCutsTable;
71}
72
[1196]73/////////////////////////////////////////////////////////////
[819]74G4ProductionCutsTable::G4ProductionCutsTable()
[961]75 : firstUse(true),verboseLevel(1),fMessenger(0)
[819]76{
77 for(size_t i=0;i< NumberOfG4CutIndex;i++)
78 {
79 rangeCutTable.push_back(new G4CutVectorForAParticle);
80 energyCutTable.push_back(new G4CutVectorForAParticle);
81 rangeDoubleVector[i] = 0;
82 energyDoubleVector[i] = 0;
83 converters[i] = 0;
84 }
85 fG4RegionStore = G4RegionStore::GetInstance();
86 defaultProductionCuts = new G4ProductionCuts();
[961]87
88 // add messenger for UI
89 fMessenger = new G4ProductionCutsTableMessenger(this);
[819]90}
91
[1196]92/////////////////////////////////////////////////////////////
[819]93G4ProductionCutsTable::G4ProductionCutsTable(const G4ProductionCutsTable& )
94{;}
95
[1196]96/////////////////////////////////////////////////////////////
[819]97G4ProductionCutsTable::~G4ProductionCutsTable()
98{
99 if (defaultProductionCuts !=0) {
100 delete defaultProductionCuts;
101 defaultProductionCuts =0;
102 }
103
104 for(CoupleTableIterator itr=coupleTable.begin();itr!=coupleTable.end();itr++){
105 delete (*itr);
106 }
107 coupleTable.clear();
108
109 for(size_t i=0;i< NumberOfG4CutIndex;i++){
110 delete rangeCutTable[i];
111 delete energyCutTable[i];
112 delete converters[i];
113 if(rangeDoubleVector[i]!=0) delete [] rangeDoubleVector[i];
114 if(energyDoubleVector[i]!=0) delete [] energyDoubleVector[i];
115 }
116 fG4ProductionCutsTable =0;
[961]117
118 if (fMessenger !=0) delete fMessenger;
119 fMessenger = 0;
[819]120}
121
122void G4ProductionCutsTable::UpdateCoupleTable(G4VPhysicalVolume* currentWorld)
123{
124 if(firstUse){
[961]125 if(G4ParticleTable::GetParticleTable()->FindParticle("gamma")){
126 converters[0] = new G4RToEConvForGamma();
127 converters[0]->SetVerboseLevel(GetVerboseLevel());
128 }
129 if(G4ParticleTable::GetParticleTable()->FindParticle("e-")){
130 converters[1] = new G4RToEConvForElectron();
131 converters[1]->SetVerboseLevel(GetVerboseLevel());
132 }
133 if(G4ParticleTable::GetParticleTable()->FindParticle("e+")){
134 converters[2] = new G4RToEConvForPositron();
135 converters[2]->SetVerboseLevel(GetVerboseLevel());
136 }
[1196]137 if(G4ParticleTable::GetParticleTable()->FindParticle("proton")){
138 converters[3] = new G4RToEConvForProton();
139 converters[3]->SetVerboseLevel(GetVerboseLevel());
140 }
[819]141 firstUse = false;
142 }
143
144 // Reset "used" flags of all couples
145 for(CoupleTableIterator CoupleItr=coupleTable.begin();
146 CoupleItr!=coupleTable.end();CoupleItr++)
147 {
148 (*CoupleItr)->SetUseFlag(false);
149 }
150
151 // Update Material-Cut-Couple
152 typedef std::vector<G4Region*>::iterator regionIterator;
153 for(regionIterator rItr=fG4RegionStore->begin();
154 rItr!=fG4RegionStore->end();rItr++)
155 {
156 // Material scan is to be done only for the regions appear in the
157 // current tracking world.
158 if((*rItr)->GetWorldPhysical()!=currentWorld) continue;
159
160 G4ProductionCuts* fProductionCut = (*rItr)->GetProductionCuts();
161 std::vector<G4Material*>::const_iterator mItr =
162 (*rItr)->GetMaterialIterator();
163 size_t nMaterial = (*rItr)->GetNumberOfMaterials();
164 (*rItr)->ClearMap();
165
166 for(size_t iMate=0;iMate<nMaterial;iMate++){
167 //check if this material cut couple has already been made
168 G4bool coupleAlreadyDefined = false;
169 G4MaterialCutsCouple* aCouple;
170 for(CoupleTableIterator cItr=coupleTable.begin();
171 cItr!=coupleTable.end();cItr++){
172 if( (*cItr)->GetMaterial()==(*mItr) &&
173 (*cItr)->GetProductionCuts()==fProductionCut){
174 coupleAlreadyDefined = true;
175 aCouple = *cItr;
176 break;
177 }
178 }
179
180 // If this combination is new, cleate and register a couple
181 if(!coupleAlreadyDefined){
182 aCouple = new G4MaterialCutsCouple((*mItr),fProductionCut);
183 coupleTable.push_back(aCouple);
184 aCouple->SetIndex(coupleTable.size()-1);
185 }
186
187 // Register this couple to the region
188 (*rItr)->RegisterMaterialCouplePair((*mItr),aCouple);
189
190 // Set the couple to the proper logical volumes in that region
191 aCouple->SetUseFlag();
192
193 std::vector<G4LogicalVolume*>::iterator rootLVItr
194 = (*rItr)->GetRootLogicalVolumeIterator();
195 size_t nRootLV = (*rItr)->GetNumberOfRootVolumes();
196 for(size_t iLV=0;iLV<nRootLV;iLV++) {
197 //Set the couple to the proper logical volumes in that region
198 G4LogicalVolume* aLV = *rootLVItr;
199 G4Region* aR = *rItr;
200
201 ScanAndSetCouple(aLV,aCouple,aR);
202
203 // Proceed to the next root logical volume in this region
204 rootLVItr++;
205 }
206
207 // Proceed to next material in this region
208 mItr++;
209 }
210 }
211
212 // Check if sizes of Range/Energy cuts tables are equal to the size of
213 // the couple table
214 // If new couples are made during the previous procedure, nCouple becomes
215 // larger then nTable
216 size_t nCouple = coupleTable.size();
217 size_t nTable = energyCutTable[0]->size();
218 G4bool newCoupleAppears = nCouple>nTable;
219 if(newCoupleAppears) {
220 for(size_t n=nCouple-nTable;n>0;n--) {
221 for(size_t nn=0;nn< NumberOfG4CutIndex;nn++){
222 rangeCutTable[nn]->push_back(-1.);
223 energyCutTable[nn]->push_back(-1.);
224 }
225 }
226 }
227
228 // Update RangeEnergy cuts tables
229 size_t idx = 0;
[1196]230 G4Timer timer;
231 if (verboseLevel>2) {
232 timer.Start();
233 }
[819]234 for(CoupleTableIterator cItr=coupleTable.begin();
235 cItr!=coupleTable.end();cItr++){
236 G4ProductionCuts* aCut = (*cItr)->GetProductionCuts();
237 const G4Material* aMat = (*cItr)->GetMaterial();
238 if((*cItr)->IsRecalcNeeded()) {
239 for(size_t ptcl=0;ptcl< NumberOfG4CutIndex;ptcl++){
240 G4double rCut = aCut->GetProductionCut(ptcl);
241 (*(rangeCutTable[ptcl]))[idx] = rCut;
242 // if(converters[ptcl] && (*cItr)->IsUsed())
243 if(converters[ptcl]) {
244 (*(energyCutTable[ptcl]))[idx] = converters[ptcl]->Convert(rCut,aMat);
245 } else {
246 (*(energyCutTable[ptcl]))[idx] = -1.;
247 }
248 }
249 }
250 idx++;
251 }
[1196]252 if (verboseLevel>2) {
253 timer.Stop();
254 std::cout << "G4ProductionCutsTable::UpdateCoupleTable "
255 << " elapsed time for calculation of energy cuts " << G4endl;
256 std::cout << timer <<G4endl;
257 }
[819]258
259 // resize Range/Energy cuts double vectors if new couple is made
260 if(newCoupleAppears){
261 for(size_t ix=0;ix<NumberOfG4CutIndex;ix++){
262 G4double* rangeVOld = rangeDoubleVector[ix];
263 G4double* energyVOld = energyDoubleVector[ix];
264 if(rangeVOld) delete [] rangeVOld;
265 if(energyVOld) delete [] energyVOld;
266 rangeDoubleVector[ix] = new G4double[(*(rangeCutTable[ix])).size()];
267 energyDoubleVector[ix] = new G4double[(*(energyCutTable[ix])).size()];
268 }
269 }
270
271 // Update Range/Energy cuts double vectors
272 for(size_t ix=0;ix<NumberOfG4CutIndex;ix++) {
273 for(size_t ixx=0;ixx<(*(rangeCutTable[ix])).size();ixx++) {
274 rangeDoubleVector[ix][ixx] = (*(rangeCutTable[ix]))[ixx];
275 energyDoubleVector[ix][ixx] = (*(energyCutTable[ix]))[ixx];
276 }
277 }
278}
279
280
[1196]281/////////////////////////////////////////////////////////////
[819]282G4double G4ProductionCutsTable::ConvertRangeToEnergy(
[1196]283 const G4ParticleDefinition* particle,
284 const G4Material* material,
285 G4double range
[819]286 )
287{
288 // This method gives energy corresponding to range value
289 // check material
290 if (material ==0) return -1.0;
291
292 // check range
293 if (range <=0.0) return -1.0;
294
295 // check particle
[1196]296 G4int index = G4ProductionCuts::GetIndex(particle);
297
[819]298 if (index<0) {
299#ifdef G4VERBOSE
300 if (verboseLevel >1) {
301 G4cout << "G4ProductionCutsTable::ConvertRangeToEnergy" ;
302 G4cout << particle->GetParticleName() << " has no cut value " << G4endl;
303 }
304#endif
305 return -1.0;
306 }
307
308 return converters[index]->Convert(range, material);
309
310}
311
[1196]312/////////////////////////////////////////////////////////////
313void G4ProductionCutsTable::ResetConverters()
314{
315 for(size_t i=0;i< NumberOfG4CutIndex;i++){
316 if (converters[i]!=0) converters[i]->Reset();
317 }
318}
[819]319
320
[1196]321/////////////////////////////////////////////////////////////
[819]322void G4ProductionCutsTable::SetEnergyRange(G4double lowedge, G4double highedge)
323{
324 G4VRangeToEnergyConverter::SetEnergyRange(lowedge,highedge);
325}
326
[1196]327/////////////////////////////////////////////////////////////
[819]328G4double G4ProductionCutsTable::GetLowEdgeEnergy() const
329{
330 return G4VRangeToEnergyConverter::GetLowEdgeEnergy();
331}
332
[1196]333/////////////////////////////////////////////////////////////
[819]334G4double G4ProductionCutsTable::GetHighEdgeEnergy() const
335{
336 return G4VRangeToEnergyConverter::GetHighEdgeEnergy();
337}
338
339
[1196]340/////////////////////////////////////////////////////////////
[819]341void G4ProductionCutsTable::ScanAndSetCouple(G4LogicalVolume* aLV,
342 G4MaterialCutsCouple* aCouple,
343 G4Region* aRegion)
344{
345 //Check whether or not this logical volume belongs to the same region
346 if((aRegion!=0) && aLV->GetRegion()!=aRegion) return;
347
348 //Check if this particular volume has a material matched to the couple
349 if(aLV->GetMaterial()==aCouple->GetMaterial()) {
350 aLV->SetMaterialCutsCouple(aCouple);
351 }
352
353 size_t noDaughters = aLV->GetNoDaughters();
354 if(noDaughters==0) return;
355
356 //Loop over daughters with same region
357 for(size_t i=0;i<noDaughters;i++){
358 G4LogicalVolume* daughterLVol = aLV->GetDaughter(i)->GetLogicalVolume();
359 ScanAndSetCouple(daughterLVol,aCouple,aRegion);
360 }
361}
362
[1196]363/////////////////////////////////////////////////////////////
[819]364void G4ProductionCutsTable::DumpCouples() const
365{
366 G4cout << G4endl;
367 G4cout << "========= Table of registered couples =============================="
368 << G4endl;
369 for(CoupleTableIterator cItr=coupleTable.begin();
370 cItr!=coupleTable.end();cItr++) {
371 G4MaterialCutsCouple* aCouple = (*cItr);
372 G4ProductionCuts* aCut = aCouple->GetProductionCuts();
373 G4cout << G4endl;
374 G4cout << "Index : " << aCouple->GetIndex()
375 << " used in the geometry : ";
376 if(aCouple->IsUsed()) G4cout << "Yes";
377 else G4cout << "No ";
378 G4cout << " recalculation needed : ";
379 if(aCouple->IsRecalcNeeded()) G4cout << "Yes";
380 else G4cout << "No ";
381 G4cout << G4endl;
382 G4cout << " Material : " << aCouple->GetMaterial()->GetName() << G4endl;
383 G4cout << " Range cuts : "
[1196]384 << " gamma " << G4BestUnit(aCut->GetProductionCut("gamma"),"Length")
385 << " e- " << G4BestUnit(aCut->GetProductionCut("e-"),"Length")
386 << " e+ " << G4BestUnit(aCut->GetProductionCut("e+"),"Length")
387 << " proton " << G4BestUnit(aCut->GetProductionCut("proton"),"Length")
[819]388 << G4endl;
389 G4cout << " Energy thresholds : " ;
390 if(aCouple->IsRecalcNeeded()) {
391 G4cout << " is not ready to print";
392 } else {
[1196]393 G4cout << " gamma " << G4BestUnit((*(energyCutTable[0]))[aCouple->GetIndex()],"Energy")
394 << " e- " << G4BestUnit((*(energyCutTable[1]))[aCouple->GetIndex()],"Energy")
395 << " e+ " << G4BestUnit((*(energyCutTable[2]))[aCouple->GetIndex()],"Energy")
396 << " proton " << G4BestUnit((*(energyCutTable[3]))[aCouple->GetIndex()],"Energy");
[819]397 }
398 G4cout << G4endl;
399
400 if(aCouple->IsUsed()) {
401 G4cout << " Region(s) which use this couple : " << G4endl;
402 typedef std::vector<G4Region*>::iterator regionIterator;
403 for(regionIterator rItr=fG4RegionStore->begin();
404 rItr!=fG4RegionStore->end();rItr++) {
405 if (IsCoupleUsedInTheRegion(aCouple, *rItr) ){
406 G4cout << " " << (*rItr)->GetName() << G4endl;
407 }
408 }
409 }
410 }
411 G4cout << G4endl;
412 G4cout << "====================================================================" << G4endl;
413 G4cout << G4endl;
414}
415
416
[1196]417/////////////////////////////////////////////////////////////
[819]418// Store cuts and material information in files under the specified directory.
419G4bool G4ProductionCutsTable::StoreCutsTable(const G4String& dir,
420 G4bool ascii)
421{
422 if (!StoreMaterialInfo(dir, ascii)) return false;
423 if (!StoreMaterialCutsCoupleInfo(dir, ascii)) return false;
424 if (!StoreCutsInfo(dir, ascii)) return false;
425
426#ifdef G4VERBOSE
[961]427 if (verboseLevel >2) {
[819]428 G4cout << "G4ProductionCutsTable::StoreCutsTable " ;
429 G4cout << " Material/Cuts information have been succesfully stored ";
430 if (ascii) {
431 G4cout << " in Ascii mode ";
432 }else{
433 G4cout << " in Binary mode ";
434 }
435 G4cout << " under " << dir << G4endl;
436 }
437#endif
438 return true;
439}
440
[1196]441/////////////////////////////////////////////////////////////
[819]442G4bool G4ProductionCutsTable::RetrieveCutsTable(const G4String& dir,
443 G4bool ascii)
444{
445 if (!CheckForRetrieveCutsTable(dir, ascii)) return false;
446 if (!RetrieveCutsInfo(dir, ascii)) return false;
447#ifdef G4VERBOSE
[961]448 if (verboseLevel >2) {
[819]449 G4cout << "G4ProductionCutsTable::RetrieveCutsTable " ;
450 G4cout << " Material/Cuts information have been succesfully retreived ";
451 if (ascii) {
452 G4cout << " in Ascii mode ";
453 }else{
454 G4cout << " in Binary mode ";
455 }
456 G4cout << " under " << dir << G4endl;
457 }
458#endif
459 return true;
460}
461
[1196]462/////////////////////////////////////////////////////////////
[819]463// check stored material and cut values are consistent
464// with the current detector setup.
465//
466G4bool
467 G4ProductionCutsTable::CheckForRetrieveCutsTable(const G4String& directory,
468 G4bool ascii)
469{
470 G4cerr << "G4ProductionCutsTable::CheckForRetrieveCutsTable!!"<< G4endl;
471 // isNeedForRestoreCoupleInfo = false;
472 if (!CheckMaterialInfo(directory, ascii)) return false;
[961]473 if (verboseLevel >2) {
[819]474 G4cerr << "G4ProductionCutsTable::CheckMaterialInfo passed !!"<< G4endl;
475 }
476 if (!CheckMaterialCutsCoupleInfo(directory, ascii)) return false;
[961]477 if (verboseLevel >2) {
[819]478 G4cerr << "G4ProductionCutsTable::CheckMaterialCutsCoupleInfo passed !!"<< G4endl;
479 }
480 return true;
481}
482
[1196]483/////////////////////////////////////////////////////////////
[819]484// Store material information in files under the specified directory.
485//
486G4bool G4ProductionCutsTable::StoreMaterialInfo(const G4String& directory,
487 G4bool ascii)
488{
489 const G4String fileName = directory + "/" + "material.dat";
[1196]490 const G4String key = "MATERIAL-V3.0";
[819]491 std::ofstream fOut;
492
493 // open output file //
494 if (!ascii )
495 fOut.open(fileName,std::ios::out|std::ios::binary);
496 else
497 fOut.open(fileName,std::ios::out);
498
499
500 // check if the file has been opened successfully
501 if (!fOut) {
502#ifdef G4VERBOSE
503 if (verboseLevel>0) {
504 G4cerr << "G4ProductionCutsTable::StoreMaterialInfo ";
505 G4cerr << " Can not open file " << fileName << G4endl;
506 }
507#endif
508 return false;
509 }
510
511 const G4MaterialTable* matTable = G4Material::GetMaterialTable();
512 // number of materials in the table
513 G4int numberOfMaterial = matTable->size();
514
515 if (ascii) {
516 /////////////// ASCII mode /////////////////
517 // key word
518 fOut << key << G4endl;
519
520 // number of materials in the table
521 fOut << numberOfMaterial << G4endl;
522
523 fOut.setf(std::ios::scientific);
524
525 // material name and density
526 for (size_t idx=0; G4int(idx)<numberOfMaterial; ++idx){
527 fOut << std::setw(FixedStringLengthForStore)
528 << ((*matTable)[idx])->GetName();
529 fOut << std::setw(FixedStringLengthForStore)
530 << ((*matTable)[idx])->GetDensity()/(g/cm3) << G4endl;
531 }
532
533 fOut.unsetf(std::ios::scientific);
534
535 } else {
536 /////////////// Binary mode /////////////////
537 char temp[FixedStringLengthForStore];
538 size_t i;
539
540 // key word
541 for (i=0; i<FixedStringLengthForStore; ++i){
542 temp[i] = '\0';
543 }
544 for (i=0; i<key.length() && i<FixedStringLengthForStore-1; ++i){
545 temp[i]=key[i];
546 }
547 fOut.write(temp, FixedStringLengthForStore);
548
549 // number of materials in the table
550 fOut.write( (char*)(&numberOfMaterial), sizeof (G4int));
551
552 // material name and density
553 for (size_t imat=0; G4int(imat)<numberOfMaterial; ++imat){
554 G4String name = ((*matTable)[imat])->GetName();
555 G4double density = ((*matTable)[imat])->GetDensity();
556 for (i=0; i<FixedStringLengthForStore; ++i) temp[i] = '\0';
557 for (i=0; i<name.length() && i<FixedStringLengthForStore-1; ++i)
558 temp[i]=name[i];
559 fOut.write(temp, FixedStringLengthForStore);
560 fOut.write( (char*)(&density), sizeof (G4double));
561 }
562 }
563
564 fOut.close();
565 return true;
566}
567
[1196]568/////////////////////////////////////////////////////////////
[819]569// check stored material is consistent with the current detector setup.
570G4bool G4ProductionCutsTable::CheckMaterialInfo(const G4String& directory,
571 G4bool ascii)
572{
573 const G4String fileName = directory + "/" + "material.dat";
[1196]574 const G4String key = "MATERIAL-V3.0";
[819]575 std::ifstream fIn;
576
577 // open input file //
578 if (!ascii )
579 fIn.open(fileName,std::ios::in|std::ios::binary);
580 else
581 fIn.open(fileName,std::ios::in);
582
583 // check if the file has been opened successfully
584 if (!fIn) {
585#ifdef G4VERBOSE
586 if (verboseLevel >0) {
587 G4cerr << "G4ProductionCutsTable::CheckMaterialInfo ";
588 G4cerr << " Can not open file " << fileName << G4endl;
589 }
590#endif
591 return false;
592 }
593
594 char temp[FixedStringLengthForStore];
595
596 // key word
597 G4String keyword;
598 if (ascii) {
599 fIn >> keyword;
600 } else {
601 fIn.read(temp, FixedStringLengthForStore);
602 keyword = (const char*)(temp);
603 }
604 if (key!=keyword) {
605#ifdef G4VERBOSE
606 if (verboseLevel >0) {
607 G4cout << "G4ProductionCutsTable::CheckMaterialInfo ";
608 G4cout << " Key word in " << fileName << "= " << keyword ;
609 G4cout <<"( should be "<< key << ")" <<G4endl;
610 }
611#endif
612 return false;
613 }
614
615 // number of materials in the table
616 G4int nmat;
617 if (ascii) {
618 fIn >> nmat;
619 } else {
620 fIn.read( (char*)(&nmat), sizeof (G4int));
621 }
622
623 // list of material
624 for (G4int idx=0; idx<nmat ; ++idx){
625 // check eof
626 if(fIn.eof()) {
627#ifdef G4VERBOSE
628 if (verboseLevel >0) {
629 G4cout << "G4ProductionCutsTable::CheckMaterialInfo ";
630 G4cout << " encountered End of File " ;
631 G4cout << " at " << idx+1 << "th material "<< G4endl;
632 }
633#endif
634 fIn.close();
635 return false;
636 }
637
638 // check material name and density
639 char name[FixedStringLengthForStore];
640 double density;
641 if (ascii) {
642 fIn >> name >> density;
643 density *= (g/cm3);
644
645 } else {
646 fIn.read(name, FixedStringLengthForStore);
647 fIn.read((char*)(&density), sizeof (G4double));
648 }
649 if (fIn.fail()) {
650#ifdef G4VERBOSE
651 if (verboseLevel >0) {
652 G4cout << "G4ProductionCutsTable::CheckMaterialInfo ";
653 G4cout << " Bad data format ";
654 G4cout << " at " << idx+1 << "th material "<< G4endl;
655 }
656#endif
657 fIn.close();
658 return false;
659 }
660
661 G4Material* aMaterial = G4Material::GetMaterial(name);
662 if (aMaterial ==0 ) continue;
663
[1055]664 G4double ratio = std::fabs(density/aMaterial->GetDensity() );
[819]665 if ((0.999>ratio) || (ratio>1.001) ){
666#ifdef G4VERBOSE
667 if (verboseLevel >0) {
668 G4cout << "G4ProductionCutsTable::CheckMaterialInfo ";
669 G4cout << " Inconsistent material density" << G4endl;;
670 G4cout << " at " << idx+1 << "th material "<< G4endl;
671 G4cout << "Name: " << name << G4endl;
672 G4cout << "Density:" << std::setiosflags(std::ios::scientific) << density / (g/cm3) ;
673 G4cout << "(should be " << aMaterial->GetDensity()/(g/cm3)<< ")" << " [g/cm3]"<< G4endl;
674 G4cout << std::resetiosflags(std::ios::scientific);
675 }
676#endif
677 fIn.close();
678 return false;
679 }
680
681 }
682
683 fIn.close();
684 return true;
685
686}
687
688
[1196]689/////////////////////////////////////////////////////////////
[819]690// Store materialCutsCouple information in files under the specified directory.
691//
692G4bool
693 G4ProductionCutsTable::StoreMaterialCutsCoupleInfo(const G4String& directory,
694 G4bool ascii)
695{
696 const G4String fileName = directory + "/" + "couple.dat";
[1196]697 const G4String key = "COUPLE-V3.0";
[819]698 std::ofstream fOut;
699 char temp[FixedStringLengthForStore];
700
701 // open output file //
702 if (!ascii )
703 fOut.open(fileName,std::ios::out|std::ios::binary);
704 else
705 fOut.open(fileName,std::ios::out);
706
707
708 // check if the file has been opened successfully
709 if (!fOut) {
710#ifdef G4VERBOSE
711 if (verboseLevel >0) {
712 G4cerr << "G4ProductionCutsTable::StoreMaterialCutsCoupleInfo ";
713 G4cerr << " Can not open file " << fileName << G4endl;
714 }
715#endif
716 return false;
717 }
718 G4int numberOfCouples = coupleTable.size();
719 if (ascii) {
720 /////////////// ASCII mode /////////////////
721 // key word
722 fOut << std::setw(FixedStringLengthForStore) << key << G4endl;
723
724 // number of couples in the table
725 fOut << numberOfCouples << G4endl;
726 } else {
727 /////////////// Binary mode /////////////////
728 // key word
729 size_t i;
730 for (i=0; i<FixedStringLengthForStore; ++i) temp[i] = '\0';
731 for (i=0; i<key.length() && i<FixedStringLengthForStore-1; ++i)
732 temp[i]=key[i];
733 fOut.write(temp, FixedStringLengthForStore);
734
735 // number of couples in the table
736 fOut.write( (char*)(&numberOfCouples), sizeof (G4int));
737 }
738
739
740 // Loop over all couples
741 CoupleTableIterator cItr;
742 for (cItr=coupleTable.begin();cItr!=coupleTable.end();cItr++){
743 G4MaterialCutsCouple* aCouple = (*cItr);
744 G4int index = aCouple->GetIndex();
745 // cut value
746 G4ProductionCuts* aCut = aCouple->GetProductionCuts();
747 G4double cutValues[NumberOfG4CutIndex];
748 for (size_t idx=0; idx <NumberOfG4CutIndex; idx++) {
749 cutValues[idx] = aCut->GetProductionCut(idx);
750 }
751 // material/region info
752 G4String materialName = aCouple->GetMaterial()->GetName();
753 G4String regionName = "NONE";
754 if (aCouple->IsUsed()){
755 typedef std::vector<G4Region*>::iterator regionIterator;
756 for(regionIterator rItr=fG4RegionStore->begin();
757 rItr!=fG4RegionStore->end();rItr++){
758 if (IsCoupleUsedInTheRegion(aCouple, *rItr) ){
759 regionName = (*rItr)->GetName();
760 break;
761 }
762 }
763 }
764
765 if (ascii) {
766 /////////////// ASCII mode /////////////////
767 // index number
768 fOut << index << G4endl;
769
770 // material name
771 fOut << std::setw(FixedStringLengthForStore) << materialName<< G4endl;
772
773 // region name
774 fOut << std::setw(FixedStringLengthForStore) << regionName<< G4endl;
775
776 fOut.setf(std::ios::scientific);
777 // cut values
778 for (size_t idx=0; idx< NumberOfG4CutIndex; idx++) {
779 fOut << std::setw(FixedStringLengthForStore) << cutValues[idx]/(mm)
780 << G4endl;
781 }
782 fOut.unsetf(std::ios::scientific);
783
784 } else {
785 /////////////// Binary mode /////////////////
786 // index
787 fOut.write( (char*)(&index), sizeof (G4int));
788
789 // material name
790 size_t i;
791 for (i=0; i<FixedStringLengthForStore; ++i) temp[i] = '\0';
792 for (i=0; i<materialName.length() && i<FixedStringLengthForStore-1; ++i) {
793 temp[i]=materialName[i];
794 }
795 fOut.write(temp, FixedStringLengthForStore);
796
797 // region name
798 for (i=0; i<FixedStringLengthForStore; ++i) temp[i] = '\0';
799 for (i=0; i<regionName.length() && i<FixedStringLengthForStore-1; ++i) {
800 temp[i]=regionName[i];
801 }
802 fOut.write(temp, FixedStringLengthForStore);
803
804 // cut values
805 for (size_t idx=0; idx< NumberOfG4CutIndex; idx++) {
806 fOut.write( (char*)(&(cutValues[idx])), sizeof (G4double));
807 }
808 }
809 }
810
811 fOut.close();
812 return true;
813}
814
815
[1196]816/////////////////////////////////////////////////////////////
[819]817// check stored materialCutsCouple is consistent
818// with the current detector setup.
819//
820G4bool
821G4ProductionCutsTable::CheckMaterialCutsCoupleInfo(const G4String& directory,
822 G4bool ascii )
823{
824 const G4String fileName = directory + "/" + "couple.dat";
[1196]825 const G4String key = "COUPLE-V3.0";
[819]826 std::ifstream fIn;
827
828 // open input file //
829 if (!ascii )
830 fIn.open(fileName,std::ios::in|std::ios::binary);
831 else
832 fIn.open(fileName,std::ios::in);
833
834 // check if the file has been opened successfully
835 if (!fIn) {
836#ifdef G4VERBOSE
837 if (verboseLevel >0) {
838 G4cerr << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo ";
839 G4cerr << " Can not open file " << fileName << G4endl;
840 }
841#endif
842 return false;
843 }
844
845 char temp[FixedStringLengthForStore];
846
847 // key word
848 G4String keyword;
849 if (ascii) {
850 fIn >> keyword;
851 } else {
852 fIn.read(temp, FixedStringLengthForStore);
853 keyword = (const char*)(temp);
854 }
855 if (key!=keyword) {
856#ifdef G4VERBOSE
857 if (verboseLevel >0) {
858 G4cout << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo ";
859 G4cout << " Key word in " << fileName << "= " << keyword ;
860 G4cout <<"( should be "<< key << ")" <<G4endl;
861 }
862#endif
863 fIn.close();
864 return false;
865 }
866
867 // numberOfCouples
868 G4int numberOfCouples;
869 if (ascii) {
870 fIn >> numberOfCouples;
871 } else {
872 fIn.read( (char*)(&numberOfCouples), sizeof (G4int));
873 }
874
875 // Reset MCCIndexConversionTable
876 mccConversionTable.Reset(numberOfCouples);
877
878 // Read in couple information
879 for (G4int idx=0; idx<numberOfCouples; idx+=1){
880 // read in index
881 G4int index;
882 if (ascii) {
883 fIn >> index;
884 } else {
885 fIn.read( (char*)(&index), sizeof (G4int));
886 }
887 // read in index material name
888 char mat_name[FixedStringLengthForStore];
889 if (ascii) {
890 fIn >> mat_name;
891 } else {
892 fIn.read(mat_name, FixedStringLengthForStore);
893 }
894 // read in index and region name
895 char region_name[FixedStringLengthForStore];
896 if (ascii) {
897 fIn >> region_name;
898 } else {
899 fIn.read(region_name, FixedStringLengthForStore);
900 }
901 // cut value
902 G4double cutValues[NumberOfG4CutIndex];
903 for (size_t i=0; i< NumberOfG4CutIndex; i++) {
904 if (ascii) {
905 fIn >> cutValues[i];
906 cutValues[i] *= (mm);
907 } else {
908 fIn.read( (char*)(&(cutValues[i])), sizeof (G4double));
909 }
910 }
911
912 // Loop over all couples
913 CoupleTableIterator cItr;
914 G4bool fOK = false;
915 G4MaterialCutsCouple* aCouple =0;
916 for (cItr=coupleTable.begin();cItr!=coupleTable.end();cItr++){
917 aCouple = (*cItr);
918 // check material name
919 if ( mat_name != aCouple->GetMaterial()->GetName() ) continue;
920 // check cut values
921 G4ProductionCuts* aCut = aCouple->GetProductionCuts();
922 G4bool fRatio = true;
923 for (size_t j=0; j< NumberOfG4CutIndex; j++) {
924 G4double ratio = cutValues[j]/aCut->GetProductionCut(j);
925 fRatio = fRatio && (0.999<ratio) && (ratio<1.001) ;
926 }
927 if (!fRatio) continue;
928 // MCC matched
929 fOK = true;
930 mccConversionTable.SetNewIndex(index, aCouple->GetIndex());
931 break;
932 }
933
934#ifdef G4VERBOSE
935 // debug information
936 if (verboseLevel >1) {
937 if (fOK) {
938 G4String regionname(region_name);
939 G4Region* fRegion = 0;
940 if ( regionname != "NONE" ) fRegion = fG4RegionStore->GetRegion(region_name);
941 if ( (( regionname == "NONE" ) && (aCouple->IsUsed()) ) ||
942 (( regionname != "NONE" ) && (fRegion==0) ) ||
943 !IsCoupleUsedInTheRegion(aCouple, fRegion) ) {
944 G4cout << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo ";
945 G4cout << "A Couple is used differnt region in the current setup ";
946 G4cout << index << ": in " << fileName << G4endl;
947 G4cout << " material: " << mat_name ;
948 G4cout << " region: " << region_name << G4endl;
949 for (size_t ii=0; ii< NumberOfG4CutIndex; ii++) {
950 G4cout << "cut[" << ii << "]=" << cutValues[ii]/mm;
951 G4cout << " mm : ";
952 }
953 G4cout << G4endl;
954 } else if ( index != aCouple->GetIndex() ) {
955 G4cout << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo ";
956 G4cout << "Index of couples was modified "<< G4endl;
957 G4cout << aCouple->GetIndex() << ":" <<aCouple->GetMaterial()->GetName();
958 G4cout <<" is defined as " ;
959 G4cout << index << ":" << mat_name << " in " << fileName << G4endl;
960 } else {
961 G4cout << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo ";
962 G4cout << index << ":" << mat_name << " in " << fileName ;
963 G4cout << " is consistent with current setup" << G4endl;
964 }
965 }
966 }
967 if ((!fOK) && (verboseLevel >0)) {
968 G4cout << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo ";
969 G4cout << "Couples is not defined in the current detector setup ";
970 G4cout << index << ": in " << fileName << G4endl;
971 G4cout << " material: " << mat_name ;
972 G4cout << " region: " << region_name << G4endl;
973 for (size_t ii=0; ii< NumberOfG4CutIndex; ii++) {
974 G4cout << "cut[" << ii << "]=" << cutValues[ii]/mm;
975 G4cout << " mm : ";
976 }
977 G4cout << G4endl;
978 }
979#endif
980
981 }
982 fIn.close();
983 return true;
984}
985
986
[1196]987/////////////////////////////////////////////////////////////
[819]988// Store cut values information in files under the specified directory.
989//
990G4bool G4ProductionCutsTable::StoreCutsInfo(const G4String& directory,
991 G4bool ascii)
992{
993 const G4String fileName = directory + "/" + "cut.dat";
[1196]994 const G4String key = "CUT-V3.0";
[819]995 std::ofstream fOut;
996 char temp[FixedStringLengthForStore];
997
998 // open output file //
999 if (!ascii )
1000 fOut.open(fileName,std::ios::out|std::ios::binary);
1001 else
1002 fOut.open(fileName,std::ios::out);
1003
1004
1005 // check if the file has been opened successfully
1006 if (!fOut) {
1007 if(verboseLevel>0) {
1008 G4cerr << "G4ProductionCutsTable::StoreCutsInfo ";
1009 G4cerr << " Can not open file " << fileName << G4endl;
1010 }
1011 return false;
1012 }
1013
1014 G4int numberOfCouples = coupleTable.size();
1015 if (ascii) {
1016 /////////////// ASCII mode /////////////////
1017 // key word
1018 fOut << key << G4endl;
1019
1020 // number of couples in the table
1021 fOut << numberOfCouples << G4endl;
1022 } else {
1023 /////////////// Binary mode /////////////////
1024 // key word
1025 size_t i;
1026 for (i=0; i<FixedStringLengthForStore; ++i) temp[i] = '\0';
1027 for (i=0; i<key.length() && i<FixedStringLengthForStore-1; ++i)
1028 temp[i]=key[i];
1029 fOut.write(temp, FixedStringLengthForStore);
1030
1031 // number of couples in the table
1032 fOut.write( (char*)(&numberOfCouples), sizeof (G4int));
1033 }
1034
1035
1036 for (size_t idx=0; idx <NumberOfG4CutIndex; idx++) {
1037 const std::vector<G4double>* fRange = GetRangeCutsVector(idx);
1038 const std::vector<G4double>* fEnergy = GetEnergyCutsVector(idx);
1039 size_t i=0;
1040 // Loop over all couples
1041 CoupleTableIterator cItr;
1042 for (cItr=coupleTable.begin();cItr!=coupleTable.end();cItr++, i++){
1043 if (ascii) {
1044 /////////////// ASCII mode /////////////////
1045 fOut.setf(std::ios::scientific);
1046 fOut << std::setw(20) << (*fRange)[i]/mm ;
1047 fOut << std::setw(20) << (*fEnergy)[i]/keV << G4endl;
1048 fOut.unsetf(std::ios::scientific);
1049 } else {
1050 /////////////// Binary mode /////////////////
1051 G4double cut = (*fRange)[i];
1052 fOut.write((char*)(&cut), sizeof (G4double));
1053 cut = (*fEnergy)[i];
1054 fOut.write((char*)(&cut), sizeof (G4double));
1055 }
1056 }
1057 }
1058 fOut.close();
1059 return true;
1060}
1061
[1196]1062/////////////////////////////////////////////////////////////
[819]1063// Retrieve cut values information in files under the specified directory.
1064//
1065G4bool G4ProductionCutsTable::RetrieveCutsInfo(const G4String& directory,
1066 G4bool ascii)
1067{
1068 const G4String fileName = directory + "/" + "cut.dat";
[1196]1069 const G4String key = "CUT-V3.0";
[819]1070 std::ifstream fIn;
1071
1072 // open input file //
1073 if (!ascii )
1074 fIn.open(fileName,std::ios::in|std::ios::binary);
1075 else
1076 fIn.open(fileName,std::ios::in);
1077
1078 // check if the file has been opened successfully
1079 if (!fIn) {
1080 if (verboseLevel >0) {
1081 G4cerr << "G4ProductionCutTable::RetrieveCutsInfo ";
1082 G4cerr << " Can not open file " << fileName << G4endl;
1083 }
1084 return false;
1085 }
1086
1087 char temp[FixedStringLengthForStore];
1088
1089 // key word
1090 G4String keyword;
1091 if (ascii) {
1092 fIn >> keyword;
1093 } else {
1094 fIn.read(temp, FixedStringLengthForStore);
1095 keyword = (const char*)(temp);
1096 }
1097 if (key!=keyword) {
1098 if (verboseLevel >0) {
1099 G4cout << "G4ProductionCutTable::RetrieveCutsInfo ";
1100 G4cout << " Key word in " << fileName << "= " << keyword ;
1101 G4cout <<"( should be "<< key << ")" <<G4endl;
1102 }
1103 return false;
1104 }
1105
1106 // numberOfCouples
1107 G4int numberOfCouples;
1108 if (ascii) {
1109 fIn >> numberOfCouples;
1110 } else {
1111 fIn.read( (char*)(&numberOfCouples), sizeof (G4int));
1112 }
1113
1114 for (size_t idx=0; G4int(idx) <NumberOfG4CutIndex; idx++) {
1115 G4CutVectorForAParticle* fRange = rangeCutTable[idx];
1116 G4CutVectorForAParticle* fEnergy = energyCutTable[idx];
1117 fRange->clear();
1118 fEnergy->clear();
1119
1120 // Loop over all couples
1121 for (size_t i=0; G4int(i)< numberOfCouples; i++){
1122 G4double rcut, ecut;
1123 if (ascii) {
1124 fIn >> rcut >> ecut;
1125 rcut *= mm;
1126 ecut *= keV;
1127 } else {
1128 fIn.read((char*)(&rcut), sizeof (G4double));
1129 fIn.read((char*)(&ecut), sizeof (G4double));
1130 }
1131 if (!mccConversionTable.IsUsed(i)) continue;
1132 size_t new_index = mccConversionTable.GetIndex(i);
1133 (*fRange)[new_index] = rcut;
1134 (*fEnergy)[new_index] = ecut;
1135 }
1136 }
1137 return true;
1138}
[961]1139
[1196]1140/////////////////////////////////////////////////////////////
[961]1141// Set Verbose Level
1142// set same verbosity to all registered RangeToEnergyConverters
1143 void G4ProductionCutsTable::SetVerboseLevel(G4int value)
1144{
1145 verboseLevel = value;
1146 for (int ip=0; ip< NumberOfG4CutIndex; ip++) {
1147 if (converters[ip] !=0 ){
1148 converters[ip]->SetVerboseLevel(value);
1149 }
1150 }
1151}
1152
[1196]1153/////////////////////////////////////////////////////////////
1154G4double G4ProductionCutsTable::GetMaxEnergyCut()
1155{
1156 return G4VRangeToEnergyConverter::GetMaxEnergyCut();
1157}
1158
1159
1160/////////////////////////////////////////////////////////////
1161void G4ProductionCutsTable::SetMaxEnergyCut(G4double value)
1162{
1163 G4VRangeToEnergyConverter::SetMaxEnergyCut(value);
1164}
Note: See TracBrowser for help on using the repository browser.