Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/cuts/src/G4ProductionCutsTable.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 // G4ProductionCutsTable class implementation
 27 //
 28 // Author: M.Asai, 5 October 2002 - First implementation
 29 // Modifications: H.Kurashige, 2004-2008
 30 // --------------------------------------------------------------------
 31 
 32 #include "G4ProductionCutsTable.hh"
 33 #include "G4ProductionCuts.hh"
 34 #include "G4MCCIndexConversionTable.hh"
 35 #include "G4ProductionCutsTableMessenger.hh"
 36 #include "G4ParticleDefinition.hh"
 37 #include "G4ParticleTable.hh"
 38 #include "G4RegionStore.hh"
 39 #include "G4LogicalVolume.hh"
 40 #include "G4VPhysicalVolume.hh"
 41 #include "G4RToEConvForElectron.hh"
 42 #include "G4RToEConvForGamma.hh"
 43 #include "G4RToEConvForPositron.hh"
 44 #include "G4RToEConvForProton.hh"
 45 #include "G4MaterialTable.hh"
 46 #include "G4Material.hh"
 47 #include "G4UnitsTable.hh"
 48 
 49 #include "G4Timer.hh"
 50 #include "G4SystemOfUnits.hh"
 51 #include "G4ios.hh"
 52 #include <iomanip>                
 53 #include <fstream>       
 54 
 55 G4ProductionCutsTable* G4ProductionCutsTable::fProductionCutsTable = nullptr;
 56 
 57 // --------------------------------------------------------------------
 58 G4ProductionCutsTable* G4ProductionCutsTable::GetProductionCutsTable()
 59 { 
 60    static G4ProductionCutsTable theProductionCutsTable;
 61    if(fProductionCutsTable == nullptr)
 62    {
 63      fProductionCutsTable = &theProductionCutsTable;
 64    }
 65   return fProductionCutsTable;
 66 }
 67 
 68 // --------------------------------------------------------------------
 69 G4ProductionCutsTable::G4ProductionCutsTable()
 70 {
 71   for(std::size_t i=0; i< NumberOfG4CutIndex; ++i)
 72   {
 73     rangeCutTable.push_back(new std::vector<G4double>);
 74     energyCutTable.push_back(new std::vector<G4double>);
 75     rangeDoubleVector[i] = nullptr;
 76     energyDoubleVector[i] = nullptr;
 77     converters[i] = nullptr;
 78   }
 79   fG4RegionStore = G4RegionStore::GetInstance();
 80   defaultProductionCuts = new G4ProductionCuts();
 81 
 82   // add messenger for UI 
 83   fMessenger = new G4ProductionCutsTableMessenger(this);
 84 }
 85 
 86 // --------------------------------------------------------------------
 87 G4ProductionCutsTable::~G4ProductionCutsTable()
 88 {
 89   delete defaultProductionCuts;
 90   defaultProductionCuts = nullptr;
 91 
 92   for(auto itr=coupleTable.cbegin(); itr!=coupleTable.cend(); ++itr)
 93   {
 94     delete (*itr); 
 95   }
 96   coupleTable.clear();
 97 
 98   for(std::size_t i=0; i< NumberOfG4CutIndex; ++i)
 99   {
100     delete rangeCutTable[i];
101     delete energyCutTable[i];
102     delete converters[i];
103     if(rangeDoubleVector[i] != nullptr) delete [] rangeDoubleVector[i];
104     if(energyDoubleVector[i] != nullptr) delete [] energyDoubleVector[i];
105     rangeCutTable[i] = nullptr;
106     energyCutTable[i] = nullptr;
107     converters[i] = nullptr;
108     rangeDoubleVector[i] = nullptr;
109     energyDoubleVector[i] = nullptr;
110     if(i < 4)
111     {
112       delete userEnergyCuts[i];
113     }
114   }
115   fProductionCutsTable = nullptr;
116 
117   delete fMessenger;
118   fMessenger = nullptr;
119 }
120 
121 // --------------------------------------------------------------------
122 void G4ProductionCutsTable::SetEnergyCutVector(const std::vector<G4double>& cutE,
123                                                std::size_t idx)
124 {
125   if (idx >= 4) {
126     G4ExceptionDescription ed;
127     ed << "Wrong index= " << idx << "; it should be < 4";
128     G4Exception("G4ProductionCutsTable::SetEnergyCutVector()",
129     "CUTS0100", FatalException, ed);
130     return;
131   }
132   userEnergyCuts[idx] = new std::vector<G4double>(cutE);
133 }
134 
135 // --------------------------------------------------------------------
136 void G4ProductionCutsTable::CreateCoupleTables()
137 {
138   // Reset "used" flags of all couples
139   for(auto CoupleItr=coupleTable.cbegin();
140            CoupleItr!=coupleTable.cend(); ++CoupleItr) 
141   {
142     (*CoupleItr)->SetUseFlag(false); 
143   }
144 
145   // Update Material-Cut-Couple
146   for(auto rItr=fG4RegionStore->cbegin(); rItr!=fG4RegionStore->cend(); ++rItr)
147   {
148     // Material scan is to be done only for the regions appear in the 
149     // current tracking world.
150     //    if((*rItr)->GetWorldPhysical()!=currentWorld) continue;
151 
152     if( (*rItr)->IsInMassGeometry() || (*rItr)->IsInParallelGeometry() )
153     {
154       G4ProductionCuts* fProductionCut = (*rItr)->GetProductionCuts();
155       auto mItr = (*rItr)->GetMaterialIterator();
156       std::size_t nMaterial = (*rItr)->GetNumberOfMaterials();
157       (*rItr)->ClearMap();
158 
159       for(std::size_t iMate=0; iMate<nMaterial; ++iMate)
160       {
161         //check if this material cut couple has already been made
162         G4bool coupleAlreadyDefined = false;
163         G4MaterialCutsCouple* aCouple;
164         for(auto cItr=coupleTable.cbegin(); cItr!=coupleTable.cend(); ++cItr)
165         {
166           if( (*cItr)->GetMaterial()==(*mItr)
167            && (*cItr)->GetProductionCuts()==fProductionCut)
168           { 
169             coupleAlreadyDefined = true;
170             aCouple = *cItr;
171             break;
172           }
173         }
174       
175         // If this combination is new, cleate and register a couple
176         if(!coupleAlreadyDefined)
177         {
178           aCouple = new G4MaterialCutsCouple((*mItr),fProductionCut);
179           coupleTable.push_back(aCouple);
180           aCouple->SetIndex(G4int(coupleTable.size()-1));
181         }
182 
183         // Register this couple to the region
184         (*rItr)->RegisterMaterialCouplePair((*mItr),aCouple);
185 
186         // Set the couple to the proper logical volumes in that region
187         aCouple->SetUseFlag();
188 
189         auto rootLVItr = (*rItr)->GetRootLogicalVolumeIterator();
190         std::size_t nRootLV = (*rItr)->GetNumberOfRootVolumes();
191         for(std::size_t iLV=0; iLV<nRootLV; ++iLV)
192         {
193           // Set the couple to the proper logical volumes in that region
194           G4LogicalVolume* aLV = *rootLVItr;
195           G4Region* aR = *rItr;
196 
197           ScanAndSetCouple(aLV,aCouple,aR);
198 
199           // Proceed to the next root logical volume in this region
200           ++rootLVItr;
201         }
202 
203         // Proceed to next material in this region
204         ++mItr;
205       }
206     }
207   }
208 
209   // Check if sizes of Range/Energy cuts tables are equal to the size of
210   // the couple table. If new couples are made during the previous procedure,
211   // nCouple becomes larger then nTable
212 
213   std::size_t nCouple = coupleTable.size();
214   std::size_t nTable = energyCutTable[0]->size();
215   G4bool newCoupleAppears = nCouple>nTable;
216   if(newCoupleAppears)
217   {
218     for(std::size_t n=nCouple-nTable; n>0; --n)
219     {
220       for(std::size_t nn=0; nn< NumberOfG4CutIndex; ++nn)
221       {
222         rangeCutTable[nn]->push_back(-1.);
223         energyCutTable[nn]->push_back(-1.);
224       }
225     }
226   }
227 
228   // resize Range/Energy cuts double vectors if new couple is made
229   if(newCoupleAppears)
230   {
231     for(std::size_t ix=0; ix<NumberOfG4CutIndex; ++ix)
232     {
233       G4double* rangeVOld = rangeDoubleVector[ix];
234       G4double* energyVOld = energyDoubleVector[ix];
235       if(rangeVOld) delete [] rangeVOld;
236       if(energyVOld) delete [] energyVOld;
237       rangeDoubleVector[ix] = new G4double[(*(rangeCutTable[ix])).size()];
238       energyDoubleVector[ix] = new G4double[(*(energyCutTable[ix])).size()];
239     }
240   }
241 }
242 
243 // --------------------------------------------------------------------
244 void G4ProductionCutsTable::UpdateCoupleTable(G4VPhysicalVolume* /*currWorld*/)
245 {
246   CreateCoupleTables();
247 
248   if(firstUse)
249   {
250     if(G4ParticleTable::GetParticleTable()->FindParticle("gamma"))
251     {
252       converters[0] = new G4RToEConvForGamma();
253       converters[0]->SetVerboseLevel(GetVerboseLevel());
254     }
255     if(G4ParticleTable::GetParticleTable()->FindParticle("e-"))
256     {
257       converters[1] = new G4RToEConvForElectron();
258       converters[1]->SetVerboseLevel(GetVerboseLevel());
259     }
260     if(G4ParticleTable::GetParticleTable()->FindParticle("e+"))
261     {
262       converters[2] = new G4RToEConvForPositron();
263       converters[2]->SetVerboseLevel(GetVerboseLevel());
264     }
265     if(G4ParticleTable::GetParticleTable()->FindParticle("proton"))
266     {
267       converters[3] = new G4RToEConvForProton();
268       converters[3]->SetVerboseLevel(GetVerboseLevel());
269     }
270     firstUse = false;
271   }
272 
273   // Update RangeEnergy cuts tables
274   std::size_t idx = 0;
275   G4Timer timer;
276   if (verboseLevel>2)
277   {
278     timer.Start();
279   }
280   for(auto cItr=coupleTable.cbegin(); cItr!=coupleTable.cend(); ++cItr)
281   {
282     G4ProductionCuts* aCut = (*cItr)->GetProductionCuts();
283     const G4Material* aMat = (*cItr)->GetMaterial();
284     if((*cItr)->IsRecalcNeeded())
285     {
286       for(std::size_t ptcl=0; ptcl< NumberOfG4CutIndex; ++ptcl)
287       {
288         G4double rCut = aCut->GetProductionCut((G4int)ptcl);
289         (*(rangeCutTable[ptcl]))[idx] = rCut;
290 
291         if(nullptr != converters[ptcl])
292         {
293           // check user defined cut in energy
294           if(nullptr == userEnergyCuts[ptcl] || userEnergyCuts[ptcl]->size() <= idx)
295           {
296             (*(energyCutTable[ptcl]))[idx] = converters[ptcl]->Convert(rCut,aMat);
297           }
298           else
299           {
300             (*(energyCutTable[ptcl]))[idx] = (*(userEnergyCuts[ptcl]))[idx]; 
301           }
302         }
303         else
304         {
305           (*(energyCutTable[ptcl]))[idx] = -1.; 
306         }
307       }
308     }
309     ++idx;  
310   }
311   if (verboseLevel>2)
312   {
313     timer.Stop();
314     G4cout << "G4ProductionCutsTable::UpdateCoupleTable() - "
315            << "Elapsed time for calculation of energy cuts: " << G4endl;
316     G4cout << timer << G4endl;
317   }
318 
319   // Update Range/Energy cuts double vectors
320   for(std::size_t ix=0; ix<NumberOfG4CutIndex; ++ix)
321   {
322     for(std::size_t ixx=0; ixx<(*(rangeCutTable[ix])).size(); ++ixx)
323     {
324       rangeDoubleVector[ix][ixx] = (*(rangeCutTable[ix]))[ixx];
325       energyDoubleVector[ix][ixx] = (*(energyCutTable[ix]))[ixx];
326     }
327   }
328 }
329 
330 // --------------------------------------------------------------------
331 G4double G4ProductionCutsTable::ConvertRangeToEnergy(
332                                   const G4ParticleDefinition* particle,
333                                   const G4Material*           material, 
334                                   G4double                    range )
335 {
336   // This method gives energy corresponding to range value  
337 
338   // protection against premature call
339   if(firstUse)
340   {
341 #ifdef G4VERBOSE
342     if(verboseLevel>0)
343     {
344       G4ExceptionDescription ed;
345       ed << "Invoked prematurely before it is fully initialized.";
346       G4Exception("G4ProductionCutsTable::ConvertRangeToEnergy()",
347                   "CUTS0100", JustWarning, ed);
348     }
349 #endif
350     return -1.0;
351   }
352 
353   // check material
354   if (material == nullptr) return -1.0;
355 
356   // check range
357   if (range == 0.0) return 0.0;
358   if (range <0.0) return -1.0;
359 
360   // check particle
361   G4int index = G4ProductionCuts::GetIndex(particle);
362 
363   if (index<0 || converters[index] == nullptr)
364   {
365 #ifdef G4VERBOSE
366     if(verboseLevel>0)
367     {
368       G4ExceptionDescription ed;
369       ed << "Invoked ";
370       if(particle != nullptr)
371       { ed << "for particle <" << particle->GetParticleName() << ">."; }
372       else
373       { ed << "without valid particle pointer."; }
374       G4Exception("G4ProductionCutsTable::ConvertRangeToEnergy()",
375                   "CUTS0101", JustWarning, ed);
376     }
377 #endif
378     return -1.0;
379   }
380 
381   return converters[index]->Convert(range, material);
382 }
383 
384 // --------------------------------------------------------------------
385 void G4ProductionCutsTable::ResetConverters()
386 {}
387 
388 // --------------------------------------------------------------------
389 void G4ProductionCutsTable::SetEnergyRange(G4double lowedge, G4double highedge)
390 {
391   G4VRangeToEnergyConverter::SetEnergyRange(lowedge,highedge);
392 }
393 
394 // --------------------------------------------------------------------
395 G4double  G4ProductionCutsTable::GetLowEdgeEnergy() const
396 {
397   return G4VRangeToEnergyConverter::GetLowEdgeEnergy();
398 }
399 
400 // --------------------------------------------------------------------
401 G4double G4ProductionCutsTable::GetHighEdgeEnergy() const
402 {
403   return G4VRangeToEnergyConverter::GetHighEdgeEnergy();
404 }
405  
406 // --------------------------------------------------------------------
407 void G4ProductionCutsTable::ScanAndSetCouple(G4LogicalVolume* aLV,
408                                              G4MaterialCutsCouple* aCouple,
409                                              G4Region* aRegion)
410 {
411   // Check whether or not this logical volume belongs to the same region
412   if((aRegion!=nullptr) && aLV->GetRegion()!=aRegion) return;
413 
414   // Check if this particular volume has a material matched to the couple
415   if(aLV->GetMaterial()==aCouple->GetMaterial())
416   {
417     aLV->SetMaterialCutsCouple(aCouple);
418   }
419 
420   std::size_t noDaughters = aLV->GetNoDaughters();
421   if(noDaughters==0) return;
422 
423   // Loop over daughters with same region
424   for(std::size_t i=0; i<noDaughters; ++i)
425   {
426     G4LogicalVolume* daughterLVol = aLV->GetDaughter(i)->GetLogicalVolume();
427     ScanAndSetCouple(daughterLVol,aCouple,aRegion);
428   }
429 }
430 
431 // --------------------------------------------------------------------
432 void G4ProductionCutsTable::DumpCouples() const
433 {
434   G4cout << G4endl;
435   G4cout << "========= Table of registered couples ============================"
436          << G4endl;
437   for(auto cItr=coupleTable.cbegin(); cItr!=coupleTable.cend(); ++cItr)
438   {
439     G4MaterialCutsCouple* aCouple = (*cItr);
440     G4ProductionCuts* aCut = aCouple->GetProductionCuts();
441     G4cout << G4endl;
442     G4cout << "Index : " << aCouple->GetIndex() 
443            << "     used in the geometry : ";
444     if(aCouple->IsUsed()) G4cout << "Yes";
445     else                  G4cout << "No ";
446 ////    G4cout << "     recalculation needed : ";
447 ////    if(aCouple->IsRecalcNeeded()) G4cout << "Yes";
448 ////    else                          G4cout << "No ";
449     G4cout << G4endl;
450     G4cout << " Material : " << aCouple->GetMaterial()->GetName() << G4endl;
451     G4cout << " Range cuts        : " 
452            << " gamma  " << G4BestUnit(aCut->GetProductionCut("gamma"),"Length")
453            << "    e-  " << G4BestUnit(aCut->GetProductionCut("e-"),"Length")
454            << "    e+  " << G4BestUnit(aCut->GetProductionCut("e+"),"Length")
455            << " proton " << G4BestUnit(aCut->GetProductionCut("proton"),"Length")
456            << G4endl;
457     G4cout << " Energy thresholds : " ;
458 ////    if(aCouple->IsRecalcNeeded()) {
459 ////      G4cout << " is not ready to print";
460 ////    } else {
461     G4cout << " gamma  " << G4BestUnit((*(energyCutTable[0]))[aCouple->GetIndex()],"Energy")
462            << "    e-  " << G4BestUnit((*(energyCutTable[1]))[aCouple->GetIndex()],"Energy")
463            << "    e+  " << G4BestUnit((*(energyCutTable[2]))[aCouple->GetIndex()],"Energy") 
464            << " proton " << G4BestUnit((*(energyCutTable[3]))[aCouple->GetIndex()],"Energy");
465 ////    }
466     G4cout << G4endl;
467 
468     if(aCouple->IsUsed())
469     {
470       G4cout << " Region(s) which use this couple : " << G4endl;
471       for(auto rItr=fG4RegionStore->cbegin();
472                rItr!=fG4RegionStore->cend(); ++rItr)
473       {
474         if (IsCoupleUsedInTheRegion(aCouple, *rItr) )
475         {
476           G4cout << "    " << (*rItr)->GetName() << G4endl;
477         }
478       }
479     }
480   }
481   G4cout << G4endl;
482   G4cout << "==================================================================" << G4endl;
483   G4cout << G4endl;
484 }
485            
486 // --------------------------------------------------------------------
487 G4bool  G4ProductionCutsTable::StoreCutsTable(const G4String& dir, 
488                                               G4bool ascii)
489 {
490   // Store cuts and material information in files under the specified directory
491 
492   if (!StoreMaterialInfo(dir, ascii)) return false;
493   if (!StoreMaterialCutsCoupleInfo(dir, ascii)) return false;
494   if (!StoreCutsInfo(dir, ascii)) return false;
495   
496 #ifdef G4VERBOSE  
497   if (verboseLevel >2)
498   {
499     G4cout << "G4ProductionCutsTable::StoreCutsTable()" << G4endl;
500     G4cout << " Material/Cuts information have been successfully stored ";
501     if (ascii)
502     {
503       G4cout << " in Ascii mode ";
504     }
505     else
506     {
507       G4cout << " in Binary mode ";
508     }
509     G4cout << " under " << dir << G4endl;  
510   }  
511 #endif
512   return true;
513 }
514   
515 // --------------------------------------------------------------------
516 G4bool  G4ProductionCutsTable::RetrieveCutsTable(const G4String& dir,
517                                                  G4bool ascii)
518 {
519   if (!CheckForRetrieveCutsTable(dir, ascii)) return false;
520   if (!RetrieveCutsInfo(dir, ascii)) return false;
521 #ifdef G4VERBOSE  
522   if (verboseLevel >2)
523   {
524     G4cout << "G4ProductionCutsTable::RetrieveCutsTable()" << G4endl;
525     G4cout << " Material/Cuts information have been successfully retrieved ";
526     if (ascii)
527     {
528       G4cout << " in Ascii mode ";
529     }
530     else
531     {
532       G4cout << " in Binary mode ";
533     }
534     G4cout << " under " << dir << G4endl;  
535   }  
536 #endif
537   return true;
538 }
539 
540 // --------------------------------------------------------------------
541 G4bool
542 G4ProductionCutsTable::CheckForRetrieveCutsTable(const G4String& directory, 
543                                                  G4bool ascii)
544 {
545   // check stored material and cut values are consistent
546   // with the current detector setup
547 
548   G4cerr << "G4ProductionCutsTable::CheckForRetrieveCutsTable()"<< G4endl;
549   //  isNeedForRestoreCoupleInfo = false;
550   if (!CheckMaterialInfo(directory, ascii)) return false;
551   if (verboseLevel >2)
552   {
553       G4cerr << "G4ProductionCutsTable::CheckMaterialInfo passed !!"<< G4endl;
554   }
555   if (!CheckMaterialCutsCoupleInfo(directory, ascii)) return false;
556   if (verboseLevel >2)
557   {
558     G4cerr << "G4ProductionCutsTable::CheckMaterialCutsCoupleInfo passed !!"
559            << G4endl;
560   }
561   return true;
562 }
563   
564 // --------------------------------------------------------------------
565 G4bool G4ProductionCutsTable::StoreMaterialInfo(const G4String& directory, 
566                                                 G4bool ascii)
567 {
568   // Store material information in files under the specified directory
569 
570   const G4String fileName = directory + "/" + "material.dat";
571   const G4String key = "MATERIAL-V3.0";
572   std::ofstream fOut;  
573 
574   // open output file
575   if (!ascii )  fOut.open(fileName,std::ios::out|std::ios::binary);
576   else          fOut.open(fileName,std::ios::out);
577 
578   // check if the file has been opened successfully 
579   if (!fOut)
580   {
581 #ifdef G4VERBOSE
582     if (verboseLevel>0)
583     {
584       G4cerr << "G4ProductionCutsTable::StoreMaterialInfo() - ";
585       G4cerr << "Cannot open file: " << fileName << G4endl;
586     }
587 #endif
588     G4Exception( "G4ProductionCutsTable::StoreMaterialInfo()",
589                  "ProcCuts102", JustWarning, "Cannot open file!");
590     return false;
591   }
592   
593   const G4MaterialTable* matTable = G4Material::GetMaterialTable(); 
594   // number of materials in the table
595   G4int numberOfMaterial = (G4int)matTable->size();
596 
597   if (ascii)
598   {
599     /////////////// ASCII mode  /////////////////
600     // key word
601     fOut  << key << G4endl;
602     
603     // number of materials in the table
604     fOut  << numberOfMaterial << G4endl;
605     
606     fOut.setf(std::ios::scientific);
607   
608     // material name and density
609     for (std::size_t idx=0; static_cast<G4int>(idx)<numberOfMaterial; ++idx)
610     {
611       fOut << std::setw(FixedStringLengthForStore)
612            << ((*matTable)[idx])->GetName();
613       fOut << std::setw(FixedStringLengthForStore)
614            << ((*matTable)[idx])->GetDensity()/(g/cm3) << G4endl;
615     }
616     
617     fOut.unsetf(std::ios::scientific);
618 
619   }
620   else
621   {
622     /////////////// Binary mode  /////////////////
623     char temp[FixedStringLengthForStore];
624     std::size_t i;
625 
626     // key word
627     for (i=0; i<FixedStringLengthForStore; ++i)
628     {
629       temp[i] = '\0'; 
630     }
631     for (i=0; i<key.length() && i<FixedStringLengthForStore-1; ++i)
632     {
633       temp[i]=key[(G4int)i];
634     }
635     fOut.write(temp, FixedStringLengthForStore);
636 
637     // number of materials in the table
638     fOut.write( (char*)(&numberOfMaterial), sizeof(G4int));
639     
640     // material name and density
641     for (std::size_t imat=0; static_cast<G4int>(imat)<numberOfMaterial; ++imat)
642     {
643       G4String name =  ((*matTable)[imat])->GetName();
644       G4double density = ((*matTable)[imat])->GetDensity();
645       for (i=0; i<FixedStringLengthForStore; ++i)
646         temp[i] = '\0'; 
647       for (i=0; i<name.length() && i<FixedStringLengthForStore-1; ++i)
648         temp[i]=name[(G4int)i];
649       fOut.write(temp, FixedStringLengthForStore);
650       fOut.write( (char*)(&density), sizeof(G4double));
651     }    
652   }    
653 
654   fOut.close();
655   return true;
656 }
657 
658 // --------------------------------------------------------------------
659 G4bool G4ProductionCutsTable::CheckMaterialInfo(const G4String& directory, 
660                                                 G4bool ascii)
661 {
662   // Check stored material is consistent with the current detector setup
663 
664   const G4String fileName = directory + "/" + "material.dat";
665   const G4String key = "MATERIAL-V3.0";
666   std::ifstream fIn;  
667 
668   // open input file
669   if (!ascii ) fIn.open(fileName,std::ios::in|std::ios::binary);
670   else         fIn.open(fileName,std::ios::in);
671 
672   // check if the file has been opened successfully 
673   if (!fIn)
674   {
675 #ifdef G4VERBOSE
676     if (verboseLevel >0)
677     {
678       G4cerr << "G4ProductionCutsTable::CheckMaterialInfo() - ";
679       G4cerr << "Cannot open file: " << fileName << G4endl;
680     }
681 #endif
682     G4Exception( "G4ProductionCutsTable::CheckMaterialInfo()",
683                  "ProcCuts102", JustWarning, "Cannot open file!");
684     return false;
685   }
686   
687   char temp[FixedStringLengthForStore];
688 
689   // key word
690   G4String keyword;    
691   if (ascii)
692   {
693     fIn >> keyword;
694   }
695   else
696   {
697     fIn.read(temp, FixedStringLengthForStore);
698     keyword = (const char*)(temp);
699   }
700   if (key!=keyword)
701   {
702 #ifdef G4VERBOSE
703     if (verboseLevel >0)
704     {
705       G4cerr << "G4ProductionCutsTable::CheckMaterialInfo() - ";
706       G4cerr << "Key word in " << fileName << "= " << keyword ;
707       G4cerr <<"( should be   "<< key << ")" <<G4endl;
708     }
709 #endif
710     G4Exception( "G4ProductionCutsTable::CheckMaterialInfo()",
711                  "ProcCuts103", JustWarning, "Bad Data Format");
712     return false;
713   }
714 
715   // number of materials in the table
716   G4int nmat;
717   if (ascii)
718   {
719     fIn >> nmat;
720   }
721   else
722   {
723     fIn.read( (char*)(&nmat), sizeof(G4int));
724   }
725   if ((nmat<=0) || (nmat >100000))
726   {
727     G4Exception( "G4ProductionCutsTable::CheckMaterialInfo()",
728                  "ProcCuts108", JustWarning, 
729                  "Number of materials is less than zero or too big");
730     return false;
731   }
732 
733   // list of material
734   for (G4int idx=0; idx<nmat ; ++idx)
735   {
736     // check eof
737     if(fIn.eof())
738     {
739 #ifdef G4VERBOSE
740       if (verboseLevel >0)
741       {
742         G4cout << "G4ProductionCutsTable::CheckMaterialInfo() - ";
743         G4cout << "Encountered End of File " ;
744         G4cout << " at " << idx+1 << "th  material "<< G4endl;
745       }
746 #endif
747       fIn.close();
748       return false;
749     }
750 
751     // check material name and density
752     char name[FixedStringLengthForStore];
753     G4double density;
754     if (ascii)
755     {
756       fIn >> name >> density;
757       density *= (g/cm3);
758       
759     }
760     else
761     {
762       fIn.read(name, FixedStringLengthForStore);
763       fIn.read((char*)(&density), sizeof(G4double));
764     }
765     if (fIn.fail())
766     {
767 #ifdef G4VERBOSE
768       if (verboseLevel >0)
769       {
770         G4cerr << "G4ProductionCutsTable::CheckMaterialInfo() - ";
771         G4cerr << "Bad data format ";
772         G4cerr << " at " << idx+1 << "th  material "<< G4endl;        
773       }
774 #endif
775       G4Exception( "G4ProductionCutsTable::CheckMaterialInfo()",
776                    "ProcCuts103", JustWarning, "Bad Data Format");
777       fIn.close();
778       return false;
779     }
780 
781     G4Material* aMaterial = G4Material::GetMaterial(name);
782     if (aMaterial == nullptr ) continue;
783 
784     G4double ratio = std::fabs(density/aMaterial->GetDensity() );
785     if ((0.999>ratio) || (ratio>1.001) )
786     {
787 #ifdef G4VERBOSE
788       if (verboseLevel >0)
789       {
790         G4cerr << "G4ProductionCutsTable::CheckMaterialInfo() - ";
791         G4cerr << "Inconsistent material density" << G4endl;;
792         G4cerr << " at " << idx+1 << "th  material "<< G4endl;        
793         G4cerr << "Name:   " << name << G4endl;
794         G4cerr << "Density:" << std::setiosflags(std::ios::scientific)
795                << density / (g/cm3) ;
796         G4cerr << "(should be " << aMaterial->GetDensity()/(g/cm3)<< ")"
797                << " [g/cm3]"<< G4endl;      
798         G4cerr << std::resetiosflags(std::ios::scientific);
799       }
800 #endif
801       G4Exception( "G4ProductionCutsTable::CheckMaterialInfo()",
802                    "ProcCuts104", JustWarning, "Inconsistent material density");
803       fIn.close();
804       return false;
805     }
806   }
807 
808   fIn.close();
809   return true;
810 }
811   
812 // --------------------------------------------------------------------
813 G4bool
814 G4ProductionCutsTable::StoreMaterialCutsCoupleInfo(const G4String& directory, 
815                                                    G4bool ascii)
816 {  
817   // Store materialCutsCouple information in files under the specified directory
818 
819   const G4String fileName = directory + "/" + "couple.dat";
820   const G4String key = "COUPLE-V3.0";
821   std::ofstream fOut;  
822   char temp[FixedStringLengthForStore];
823 
824   // open output file
825   if (!ascii ) fOut.open(fileName,std::ios::out|std::ios::binary);
826   else         fOut.open(fileName,std::ios::out);
827   
828   
829   // check if the file has been opened successfully 
830   if (!fOut)
831   {
832 #ifdef G4VERBOSE
833     if (verboseLevel >0)
834     {
835       G4cerr << "G4ProductionCutsTable::StoreMaterialCutsCoupleInfo() - ";
836       G4cerr << "Cannot open file: " << fileName << G4endl;
837     }
838 #endif
839     G4Exception( "G4ProductionCutsTable::StoreMaterialCutsCoupleInfo()",
840                  "ProcCuts102",
841                  JustWarning, "Cannot open file!");
842      return false;
843   }
844   G4int numberOfCouples = (G4int)coupleTable.size();
845   if (ascii)
846   {
847     /////////////// ASCII mode  /////////////////
848     // key word
849     fOut << std::setw(FixedStringLengthForStore) <<  key << G4endl;
850     
851     // number of couples in the table
852     fOut  << numberOfCouples << G4endl;
853   }
854   else
855   {
856     /////////////// Binary mode  /////////////////
857     // key word
858     std::size_t i;
859     for (i=0; i<FixedStringLengthForStore; ++i)
860       temp[i] = '\0'; 
861     for (i=0; i<key.length() && i<FixedStringLengthForStore-1; ++i)
862       temp[i]=key[(G4int)i];
863     fOut.write(temp, FixedStringLengthForStore);
864     
865     // number of couples in the table   
866     fOut.write( (char*)(&numberOfCouples), sizeof(G4int));
867   }
868 
869   // Loop over all couples
870   for (auto cItr=coupleTable.cbegin(); cItr!=coupleTable.cend(); ++cItr)
871   {
872     G4MaterialCutsCouple* aCouple = (*cItr);
873     G4int index = aCouple->GetIndex(); 
874     // cut value
875     G4ProductionCuts* aCut = aCouple->GetProductionCuts();
876     G4double cutValues[NumberOfG4CutIndex];
877     for (std::size_t idx=0; idx <NumberOfG4CutIndex; ++idx)
878     {
879       cutValues[idx] = aCut->GetProductionCut((G4int)idx);
880     }
881     // material/region info
882     G4String materialName = aCouple->GetMaterial()->GetName();
883     G4String regionName = "NONE";
884     if (aCouple->IsUsed())
885     {
886       for(auto rItr=fG4RegionStore->cbegin();
887                rItr!=fG4RegionStore->cend(); ++rItr)
888       {
889         if (IsCoupleUsedInTheRegion(aCouple, *rItr))
890         {
891           regionName = (*rItr)->GetName();
892           break;
893         }
894       }
895     }
896 
897     if (ascii)
898     {
899       /////////////// ASCII mode  /////////////////
900       // index number
901       fOut  << index << G4endl; 
902   
903       // material name 
904       fOut << std::setw(FixedStringLengthForStore) << materialName<< G4endl;
905   
906       // region name 
907       fOut << std::setw(FixedStringLengthForStore) << regionName<< G4endl;
908 
909       fOut.setf(std::ios::scientific);
910       // cut values
911       for (std::size_t idx=0; idx< NumberOfG4CutIndex; ++idx)
912       {
913         fOut << std::setw(FixedStringLengthForStore) << cutValues[idx]/(mm)
914              << G4endl;
915       }
916       fOut.unsetf(std::ios::scientific);
917 
918     }
919     else
920     {
921       /////////////// Binary mode  /////////////////
922       // index
923       fOut.write( (char*)(&index), sizeof(G4int));
924     
925       // material name
926       std::size_t i;
927       for (i=0; i<FixedStringLengthForStore; ++i)
928         temp[i] = '\0'; 
929       for (i=0; i<materialName.length() && i<FixedStringLengthForStore-1; ++i)
930         temp[i]=materialName[(G4int)i];
931       fOut.write(temp, FixedStringLengthForStore);
932 
933       // region name
934       for (i=0; i<FixedStringLengthForStore; ++i)
935         temp[i] = '\0'; 
936       for (i=0; i<regionName.length() && i<FixedStringLengthForStore-1; ++i)
937         temp[i]=regionName[(G4int)i];
938       fOut.write(temp, FixedStringLengthForStore);
939 
940       // cut values
941       for (std::size_t idx=0; idx< NumberOfG4CutIndex; ++idx)
942       {
943          fOut.write( (char*)(&(cutValues[idx])), sizeof(G4double));
944       }    
945     }
946   }
947   fOut.close();
948   return true;
949 }
950 
951 // --------------------------------------------------------------------
952 G4bool
953 G4ProductionCutsTable::CheckMaterialCutsCoupleInfo(const G4String& directory,
954                                                    G4bool ascii)
955 {
956   // Check stored materialCutsCouple is consistent
957   // with the current detector setup. 
958 
959   const G4String fileName = directory + "/" + "couple.dat";
960   const G4String key = "COUPLE-V3.0";
961   std::ifstream fIn;  
962 
963   // open input file
964   if (!ascii ) fIn.open(fileName,std::ios::in|std::ios::binary);
965   else         fIn.open(fileName,std::ios::in);
966 
967   // check if the file has been opened successfully 
968   if (!fIn)
969   {
970 #ifdef G4VERBOSE
971     if (verboseLevel >0)
972     {
973       G4cerr << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo() - ";
974       G4cerr << "Cannot open file!" << fileName << G4endl;
975     }
976 #endif
977     G4Exception( "G4ProductionCutsTable::CheckMaterialCutsCoupleInfo()",
978                  "ProcCuts102", JustWarning, "Cannot open file!");
979     return false;
980   }
981   
982   char temp[FixedStringLengthForStore];
983 
984    // key word
985   G4String keyword;    
986   if (ascii)
987   {
988     fIn >> keyword;
989   }
990   else
991   {
992     fIn.read(temp, FixedStringLengthForStore);
993     keyword = (const char*)(temp);
994   }
995   if (key!=keyword)
996   {
997 #ifdef G4VERBOSE
998     if (verboseLevel >0)
999     {
1000       G4cerr << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo() - ";
1001       G4cerr << "Key word in " << fileName << "= " << keyword ;
1002       G4cerr <<"( should be   "<< key << ")" << G4endl;
1003     }
1004 #endif
1005     G4Exception( "G4ProductionCutsTable::CheckMaterialCutsCoupleInfo()",
1006                  "ProcCuts103", JustWarning, "Bad Data Format");
1007     fIn.close();
1008     return false;
1009   }
1010 
1011   // numberOfCouples
1012   G4int numberOfCouples;    
1013   if (ascii)
1014   {
1015     fIn >> numberOfCouples;
1016   }
1017   else
1018   {
1019     fIn.read( (char*)(&numberOfCouples), sizeof(G4int));
1020   }
1021 
1022   // Reset MCCIndexConversionTable
1023   mccConversionTable.Reset(numberOfCouples);  
1024 
1025   // Read in couple information
1026   for (G4int idx=0; idx<numberOfCouples; ++idx)
1027   {
1028     // read in index
1029     G4int index; 
1030     if (ascii)
1031     {
1032       fIn >> index;
1033     }
1034     else
1035     {
1036       fIn.read( (char*)(&index), sizeof(G4int));
1037     }
1038     // read in index material name
1039     char mat_name[FixedStringLengthForStore];
1040     if (ascii)
1041     {
1042       fIn >> mat_name;
1043     }
1044     else
1045     {
1046       fIn.read(mat_name, FixedStringLengthForStore);
1047     }
1048     // read in index and region name
1049     char region_name[FixedStringLengthForStore];
1050     if (ascii)
1051     {
1052       fIn >> region_name;
1053     }
1054     else
1055     {
1056       fIn.read(region_name, FixedStringLengthForStore);
1057     }
1058     // cut value
1059     G4double cutValues[NumberOfG4CutIndex];
1060     for (std::size_t i=0; i< NumberOfG4CutIndex; ++i)
1061     {
1062       if (ascii)
1063       {
1064         fIn >>  cutValues[i];
1065         cutValues[i] *= (mm);
1066       }
1067       else
1068       {
1069         fIn.read( (char*)(&(cutValues[i])), sizeof(G4double));
1070       }
1071     }
1072  
1073     // Loop over all couples
1074     G4bool fOK = false;
1075     G4MaterialCutsCouple* aCouple = nullptr;
1076     for (auto cItr=coupleTable.cbegin(); cItr!=coupleTable.cend(); ++cItr)
1077     {
1078       aCouple = (*cItr);
1079       // check material name
1080       if ( mat_name !=  aCouple->GetMaterial()->GetName() ) continue;
1081       // check cut values
1082       G4ProductionCuts* aCut = aCouple->GetProductionCuts();
1083       G4bool fRatio = true;
1084       for (std::size_t j=0; j< NumberOfG4CutIndex; ++j)
1085       {
1086         // check ratio only if values are not the same
1087         if (cutValues[j] != aCut->GetProductionCut((G4int)j))
1088         {  
1089           G4double ratio =  cutValues[j]/aCut->GetProductionCut((G4int)j);
1090           fRatio = fRatio && (0.999<ratio) && (ratio<1.001) ;
1091         }
1092       } 
1093       if (!fRatio) continue; 
1094       // MCC matched 
1095       fOK = true;
1096       mccConversionTable.SetNewIndex(index, aCouple->GetIndex()); 
1097       break;
1098     }
1099 
1100     if (fOK)
1101     {
1102 #ifdef G4VERBOSE
1103       // debug information 
1104       if (verboseLevel >1)
1105       {
1106         G4String regionname(region_name);
1107         G4Region* fRegion = nullptr;
1108         if ( regionname != "NONE" )
1109         {
1110           fRegion = fG4RegionStore->GetRegion(region_name);
1111           if (fRegion == nullptr)
1112           {
1113             G4cout << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo() - ";
1114             G4cout << "Region " << regionname << " is not found ";            
1115             G4cout << index << ": in " << fileName  << G4endl;
1116           } 
1117         }
1118         if  (((regionname == "NONE") && (aCouple->IsUsed()))
1119          || ((fRegion!=nullptr) && !IsCoupleUsedInTheRegion(aCouple, fRegion)))
1120         {
1121           G4cout << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo()"
1122                  << G4endl;
1123           G4cout << "A Couple is used different region in the current setup ";
1124           G4cout << index << ": in " << fileName  << G4endl;
1125           G4cout << " material: " << mat_name ;
1126           G4cout << " region: " << region_name << G4endl;
1127           for (std::size_t ii=0; ii< NumberOfG4CutIndex; ++ii)
1128           {
1129             G4cout << "cut[" << ii << "]=" << cutValues[ii]/mm;
1130             G4cout << " mm   :  ";
1131           } 
1132           G4cout << G4endl;
1133         }
1134         else if ( index !=  aCouple->GetIndex() )
1135         {
1136           G4cout << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo() - ";
1137           G4cout << "Index of couples was modified "<< G4endl;
1138           G4cout << aCouple->GetIndex() << ":"
1139                  << aCouple->GetMaterial()->GetName();
1140           G4cout <<" is defined as " ;
1141           G4cout << index << ":"  << mat_name << " in " << fileName << G4endl;
1142         }
1143         else
1144         {
1145           G4cout << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo() - ";
1146           G4cout << index << ":"  << mat_name << " in " << fileName ;
1147           G4cout << " is consistent with current setup" << G4endl;
1148         }
1149       }
1150 #endif
1151     }
1152     else
1153     {
1154 #ifdef G4VERBOSE
1155       if (verboseLevel >0)
1156       {
1157         G4cout << "G4ProductionCutTable::CheckMaterialCutsCoupleInfo()"
1158                << G4endl;
1159         G4cout << "Couples are not defined in the current detector setup ";
1160         G4cout << index << ": in " << fileName  << G4endl;
1161         G4cout << " material: " << mat_name ;
1162         G4cout << " region: " << region_name << G4endl;
1163         for (std::size_t ii=0; ii< NumberOfG4CutIndex; ++ii)
1164         {
1165           G4cout << "cut[" << ii << "]=" << cutValues[ii]/mm;
1166           G4cout << " mm   :  ";
1167         }
1168         G4cout << G4endl;
1169       }
1170 #endif
1171     }
1172   }
1173   fIn.close();
1174   return true;
1175 }
1176 
1177 // --------------------------------------------------------------------
1178 G4bool G4ProductionCutsTable::StoreCutsInfo(const G4String& directory, 
1179                                             G4bool ascii)
1180 {
1181   // Store cut values information in files under the specified directory
1182 
1183   const G4String fileName = directory + "/" + "cut.dat";
1184   const G4String key = "CUT-V3.0";
1185   std::ofstream fOut;  
1186   char temp[FixedStringLengthForStore];
1187   
1188   // open output file
1189   if (!ascii ) fOut.open(fileName,std::ios::out|std::ios::binary);
1190   else         fOut.open(fileName,std::ios::out);
1191 
1192   // check if the file has been opened successfully 
1193   if (!fOut)
1194   {
1195     if(verboseLevel>0)
1196     {         
1197       G4cerr << "G4ProductionCutsTable::StoreCutsInfo() - ";
1198       G4cerr << "Cannot open file: " << fileName << G4endl;
1199     }
1200     G4Exception( "G4ProductionCutsTable::StoreCutsInfo()",
1201                  "ProcCuts102", JustWarning, "Cannot open file!");
1202     return false;
1203   }
1204 
1205   G4int numberOfCouples = (G4int)coupleTable.size();
1206   if (ascii)
1207   {
1208     /////////////// ASCII mode  /////////////////
1209     // key word
1210     fOut << key << G4endl;
1211     
1212     // number of couples in the table
1213     fOut << numberOfCouples << G4endl;
1214   }
1215   else
1216   {
1217     /////////////// Binary mode  /////////////////
1218     // key word
1219     std::size_t i;
1220     for (i=0; i<FixedStringLengthForStore; ++i)
1221       temp[i] = '\0'; 
1222     for (i=0; i<key.length() && i<FixedStringLengthForStore-1; ++i)
1223       temp[i]=key[(G4int)i];
1224     fOut.write(temp, FixedStringLengthForStore);
1225     
1226     // number of couples in the table   
1227     fOut.write( (char*)(&numberOfCouples), sizeof(G4int));
1228   }
1229 
1230   for (std::size_t idx=0; idx <NumberOfG4CutIndex; ++idx)
1231   {
1232     const std::vector<G4double>* fRange  = GetRangeCutsVector(idx);
1233     const std::vector<G4double>* fEnergy = GetEnergyCutsVector(idx);
1234     std::size_t i=0;
1235     // Loop over all couples
1236     for (auto cItr=coupleTable.cbegin();cItr!=coupleTable.cend(); ++cItr, ++i)
1237     {      
1238       if (ascii)
1239       {
1240         /////////////// ASCII mode  /////////////////
1241         fOut.setf(std::ios::scientific);
1242         fOut << std::setw(20) << (*fRange)[i]/mm  ;
1243         fOut << std::setw(20) << (*fEnergy)[i]/keV << G4endl;
1244         fOut.unsetf(std::ios::scientific);
1245       }
1246       else
1247       {
1248         /////////////// Binary mode  /////////////////
1249         G4double cut =  (*fRange)[i];
1250         fOut.write((char*)(&cut), sizeof(G4double));
1251         cut =  (*fEnergy)[i];
1252         fOut.write((char*)(&cut), sizeof(G4double));
1253       }
1254     }
1255   }
1256   fOut.close();
1257   return true;
1258 }
1259   
1260 // --------------------------------------------------------------------
1261 G4bool G4ProductionCutsTable::RetrieveCutsInfo(const G4String& directory,
1262                                                G4bool ascii)
1263 {
1264   // Retrieve cut values information in files under the specified directory
1265 
1266   const G4String fileName = directory + "/" + "cut.dat";
1267   const G4String key = "CUT-V3.0";
1268   std::ifstream fIn;  
1269 
1270   // open input file
1271   if (!ascii ) fIn.open(fileName,std::ios::in|std::ios::binary);
1272   else         fIn.open(fileName,std::ios::in);
1273 
1274   // check if the file has been opened successfully 
1275   if (!fIn)
1276   {
1277     if (verboseLevel >0)
1278     {
1279       G4cerr << "G4ProductionCutTable::RetrieveCutsInfo() - ";
1280       G4cerr << "Cannot open file: " << fileName << G4endl;
1281     }
1282     G4Exception( "G4ProductionCutsTable::RetrieveCutsInfo()",
1283                  "ProcCuts102", JustWarning, "Cannot open file!");
1284     return false;
1285   }
1286   
1287   char temp[FixedStringLengthForStore];
1288 
1289    // key word
1290   G4String keyword;    
1291   if (ascii)
1292   {
1293     fIn >> keyword;
1294   }
1295   else
1296   {
1297     fIn.read(temp, FixedStringLengthForStore);
1298     keyword = (const char*)(temp);
1299   }
1300   if (key!=keyword)
1301   {
1302     if (verboseLevel >0)
1303     {
1304       G4cerr << "G4ProductionCutTable::RetrieveCutsInfo() - ";
1305       G4cerr << "Key word in " << fileName << "= " << keyword ;
1306       G4cerr <<"( should be   "<< key << ")" << G4endl;
1307     }
1308     G4Exception( "G4ProductionCutsTable::RetrieveCutsInfo()",
1309                  "ProcCuts103", JustWarning, "Bad Data Format");
1310     return false;
1311   }
1312 
1313   // numberOfCouples
1314   G4int numberOfCouples;    
1315   if (ascii)
1316   {
1317     fIn >> numberOfCouples;
1318     if (fIn.fail())
1319     {
1320       G4Exception( "G4ProductionCutsTable::RetrieveCutsInfo()",
1321                    "ProcCuts103", JustWarning, "Bad Data Format");
1322       return false;
1323     }
1324   }
1325   else
1326   {
1327     fIn.read( (char*)(&numberOfCouples), sizeof(G4int));
1328   }
1329 
1330   if (numberOfCouples > static_cast<G4int>(mccConversionTable.size()) )
1331   {
1332     G4Exception( "G4ProductionCutsTable::RetrieveCutsInfo()",
1333                  "ProcCuts109", JustWarning, 
1334                  "Number of Couples in the file exceeds defined couples");
1335   }
1336   numberOfCouples = (G4int)mccConversionTable.size();
1337   
1338   for (std::size_t idx=0; static_cast<G4int>(idx) <NumberOfG4CutIndex; ++idx)
1339   {
1340     std::vector<G4double>* fRange  = rangeCutTable[idx];
1341     std::vector<G4double>* fEnergy = energyCutTable[idx];
1342     fRange->clear();
1343     fEnergy->clear();
1344 
1345     // Loop over all couples
1346     for (std::size_t i=0; static_cast<G4int>(i)< numberOfCouples; ++i)
1347     {      
1348       G4double rcut, ecut;
1349       if (ascii)
1350       {
1351         fIn >> rcut >> ecut;
1352         if (fIn.fail())
1353         {
1354           G4Exception( "G4ProductionCutsTable::RetrieveCutsInfo()",
1355                        "ProcCuts103", JustWarning, "Bad Data Format");
1356           return false;
1357         }
1358         rcut *= mm;
1359         ecut *= keV;
1360       }
1361       else
1362       {
1363         fIn.read((char*)(&rcut), sizeof(G4double));
1364         fIn.read((char*)(&ecut), sizeof(G4double));
1365       }
1366       if (!mccConversionTable.IsUsed(i)) continue;
1367       std::size_t new_index = mccConversionTable.GetIndex(i);
1368       (*fRange)[new_index]  = rcut;
1369       (*fEnergy)[new_index] = ecut;
1370     }
1371   }
1372   return true;
1373 }
1374 
1375 // --------------------------------------------------------------------
1376 void G4ProductionCutsTable::SetVerboseLevel(G4int value)
1377 {
1378   // Set same verbosity to all registered RangeToEnergyConverters  
1379 
1380   verboseLevel = value;
1381   for (G4int ip=0; ip< NumberOfG4CutIndex; ++ip)
1382   {
1383     if (converters[ip] != nullptr )
1384     {
1385       converters[ip]->SetVerboseLevel(value);
1386     }
1387   }
1388 }
1389 
1390 // --------------------------------------------------------------------
1391 G4double G4ProductionCutsTable::GetMaxEnergyCut()
1392 {
1393   return G4VRangeToEnergyConverter::GetMaxEnergyCut();
1394 }
1395 
1396 // --------------------------------------------------------------------
1397 void G4ProductionCutsTable::SetMaxEnergyCut(G4double value)
1398 {
1399   G4VRangeToEnergyConverter::SetMaxEnergyCut(value);
1400 }
1401