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

Last change on this file since 833 was 819, checked in by garnier, 17 years ago

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