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

Last change on this file since 1228 was 1228, checked in by garnier, 14 years ago

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