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

Last change on this file since 1315 was 1315, checked in by garnier, 15 years ago

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