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

Last change on this file since 1315 was 1315, checked in by garnier, 14 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.