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

Last change on this file since 1138 was 1055, checked in by garnier, 17 years ago

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