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

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

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