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

Last change on this file since 869 was 819, checked in by garnier, 16 years ago

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