| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4ProductionCutsTable class implementation
27 //
28 // Author: M.Asai, 5 October 2002 - First impl
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::
56
57 // -------------------------------------------
58 G4ProductionCutsTable* G4ProductionCutsTable::
59 {
60 static G4ProductionCutsTable theProductionC
61 if(fProductionCutsTable == nullptr)
62 {
63 fProductionCutsTable = &theProductionCuts
64 }
65 return fProductionCutsTable;
66 }
67
68 // -------------------------------------------
69 G4ProductionCutsTable::G4ProductionCutsTable()
70 {
71 for(std::size_t i=0; i< NumberOfG4CutIndex;
72 {
73 rangeCutTable.push_back(new std::vector<G4
74 energyCutTable.push_back(new std::vector<G
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 G4ProductionCutsTableMessen
84 }
85
86 // -------------------------------------------
87 G4ProductionCutsTable::~G4ProductionCutsTable(
88 {
89 delete defaultProductionCuts;
90 defaultProductionCuts = nullptr;
91
92 for(auto itr=coupleTable.cbegin(); itr!=coup
93 {
94 delete (*itr);
95 }
96 coupleTable.clear();
97
98 for(std::size_t i=0; i< NumberOfG4CutIndex;
99 {
100 delete rangeCutTable[i];
101 delete energyCutTable[i];
102 delete converters[i];
103 if(rangeDoubleVector[i] != nullptr) delete
104 if(energyDoubleVector[i] != nullptr) delet
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
123
124 {
125 if (idx >= 4) {
126 G4ExceptionDescription ed;
127 ed << "Wrong index= " << idx << "; it shou
128 G4Exception("G4ProductionCutsTable::SetEne
129 "CUTS0100", FatalException, ed);
130 return;
131 }
132 userEnergyCuts[idx] = new std::vector<G4doub
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(); ++Co
141 {
142 (*CoupleItr)->SetUseFlag(false);
143 }
144
145 // Update Material-Cut-Couple
146 for(auto rItr=fG4RegionStore->cbegin(); rItr
147 {
148 // Material scan is to be done only for th
149 // current tracking world.
150 // if((*rItr)->GetWorldPhysical()!=curr
151
152 if( (*rItr)->IsInMassGeometry() || (*rItr)
153 {
154 G4ProductionCuts* fProductionCut = (*rIt
155 auto mItr = (*rItr)->GetMaterialIterator
156 std::size_t nMaterial = (*rItr)->GetNumb
157 (*rItr)->ClearMap();
158
159 for(std::size_t iMate=0; iMate<nMaterial
160 {
161 //check if this material cut couple ha
162 G4bool coupleAlreadyDefined = false;
163 G4MaterialCutsCouple* aCouple;
164 for(auto cItr=coupleTable.cbegin(); cI
165 {
166 if( (*cItr)->GetMaterial()==(*mItr)
167 && (*cItr)->GetProductionCuts()==fP
168 {
169 coupleAlreadyDefined = true;
170 aCouple = *cItr;
171 break;
172 }
173 }
174
175 // If this combination is new, cleate
176 if(!coupleAlreadyDefined)
177 {
178 aCouple = new G4MaterialCutsCouple((
179 coupleTable.push_back(aCouple);
180 aCouple->SetIndex(G4int(coupleTable.
181 }
182
183 // Register this couple to the region
184 (*rItr)->RegisterMaterialCouplePair((*
185
186 // Set the couple to the proper logica
187 aCouple->SetUseFlag();
188
189 auto rootLVItr = (*rItr)->GetRootLogic
190 std::size_t nRootLV = (*rItr)->GetNumb
191 for(std::size_t iLV=0; iLV<nRootLV; ++
192 {
193 // Set the couple to the proper logi
194 G4LogicalVolume* aLV = *rootLVItr;
195 G4Region* aR = *rItr;
196
197 ScanAndSetCouple(aLV,aCouple,aR);
198
199 // Proceed to the next root logical
200 ++rootLVItr;
201 }
202
203 // Proceed to next material in this re
204 ++mItr;
205 }
206 }
207 }
208
209 // Check if sizes of Range/Energy cuts table
210 // the couple table. If new couples are made
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< NumberOfG4CutI
221 {
222 rangeCutTable[nn]->push_back(-1.);
223 energyCutTable[nn]->push_back(-1.);
224 }
225 }
226 }
227
228 // resize Range/Energy cuts double vectors i
229 if(newCoupleAppears)
230 {
231 for(std::size_t ix=0; ix<NumberOfG4CutInde
232 {
233 G4double* rangeVOld = rangeDoubleVector[
234 G4double* energyVOld = energyDoubleVecto
235 if(rangeVOld) delete [] rangeVOld;
236 if(energyVOld) delete [] energyVOld;
237 rangeDoubleVector[ix] = new G4double[(*(
238 energyDoubleVector[ix] = new G4double[(*
239 }
240 }
241 }
242
243 // -------------------------------------------
244 void G4ProductionCutsTable::UpdateCoupleTable(
245 {
246 CreateCoupleTables();
247
248 if(firstUse)
249 {
250 if(G4ParticleTable::GetParticleTable()->Fi
251 {
252 converters[0] = new G4RToEConvForGamma()
253 converters[0]->SetVerboseLevel(GetVerbos
254 }
255 if(G4ParticleTable::GetParticleTable()->Fi
256 {
257 converters[1] = new G4RToEConvForElectro
258 converters[1]->SetVerboseLevel(GetVerbos
259 }
260 if(G4ParticleTable::GetParticleTable()->Fi
261 {
262 converters[2] = new G4RToEConvForPositro
263 converters[2]->SetVerboseLevel(GetVerbos
264 }
265 if(G4ParticleTable::GetParticleTable()->Fi
266 {
267 converters[3] = new G4RToEConvForProton(
268 converters[3]->SetVerboseLevel(GetVerbos
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!=co
281 {
282 G4ProductionCuts* aCut = (*cItr)->GetProdu
283 const G4Material* aMat = (*cItr)->GetMater
284 if((*cItr)->IsRecalcNeeded())
285 {
286 for(std::size_t ptcl=0; ptcl< NumberOfG4
287 {
288 G4double rCut = aCut->GetProductionCut
289 (*(rangeCutTable[ptcl]))[idx] = rCut;
290
291 if(nullptr != converters[ptcl])
292 {
293 // check user defined cut in energy
294 if(nullptr == userEnergyCuts[ptcl] |
295 {
296 (*(energyCutTable[ptcl]))[idx] = c
297 }
298 else
299 {
300 (*(energyCutTable[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::UpdateCo
315 << "Elapsed time for calculation of
316 G4cout << timer << G4endl;
317 }
318
319 // Update Range/Energy cuts double vectors
320 for(std::size_t ix=0; ix<NumberOfG4CutIndex;
321 {
322 for(std::size_t ixx=0; ixx<(*(rangeCutTabl
323 {
324 rangeDoubleVector[ix][ixx] = (*(rangeCut
325 energyDoubleVector[ix][ixx] = (*(energyC
326 }
327 }
328 }
329
330 // -------------------------------------------
331 G4double G4ProductionCutsTable::ConvertRangeTo
332 const G4Part
333 const G4Mate
334 G4double
335 {
336 // This method gives energy corresponding to
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
346 G4Exception("G4ProductionCutsTable::Conv
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(par
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->Ge
372 else
373 { ed << "without valid particle pointer.
374 G4Exception("G4ProductionCutsTable::Conv
375 "CUTS0101", JustWarning, ed)
376 }
377 #endif
378 return -1.0;
379 }
380
381 return converters[index]->Convert(range, mat
382 }
383
384 // -------------------------------------------
385 void G4ProductionCutsTable::ResetConverters()
386 {}
387
388 // -------------------------------------------
389 void G4ProductionCutsTable::SetEnergyRange(G4d
390 {
391 G4VRangeToEnergyConverter::SetEnergyRange(lo
392 }
393
394 // -------------------------------------------
395 G4double G4ProductionCutsTable::GetLowEdgeEne
396 {
397 return G4VRangeToEnergyConverter::GetLowEdge
398 }
399
400 // -------------------------------------------
401 G4double G4ProductionCutsTable::GetHighEdgeEne
402 {
403 return G4VRangeToEnergyConverter::GetHighEdg
404 }
405
406 // -------------------------------------------
407 void G4ProductionCutsTable::ScanAndSetCouple(G
408 G
409 G
410 {
411 // Check whether or not this logical volume
412 if((aRegion!=nullptr) && aLV->GetRegion()!=a
413
414 // Check if this particular volume has a mat
415 if(aLV->GetMaterial()==aCouple->GetMaterial(
416 {
417 aLV->SetMaterialCutsCouple(aCouple);
418 }
419
420 std::size_t noDaughters = aLV->GetNoDaughter
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->GetDa
427 ScanAndSetCouple(daughterLVol,aCouple,aReg
428 }
429 }
430
431 // -------------------------------------------
432 void G4ProductionCutsTable::DumpCouples() cons
433 {
434 G4cout << G4endl;
435 G4cout << "========= Table of registered cou
436 << G4endl;
437 for(auto cItr=coupleTable.cbegin(); cItr!=co
438 {
439 G4MaterialCutsCouple* aCouple = (*cItr);
440 G4ProductionCuts* aCut = aCouple->GetProdu
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 <
448 //// else G4cout <
449 G4cout << G4endl;
450 G4cout << " Material : " << aCouple->GetMa
451 G4cout << " Range cuts : "
452 << " gamma " << G4BestUnit(aCut->G
453 << " e- " << G4BestUnit(aCut->G
454 << " e+ " << G4BestUnit(aCut->G
455 << " proton " << G4BestUnit(aCut->G
456 << G4endl;
457 G4cout << " Energy thresholds : " ;
458 //// if(aCouple->IsRecalcNeeded()) {
459 //// G4cout << " is not ready to print";
460 //// } else {
461 G4cout << " gamma " << G4BestUnit((*(ener
462 << " e- " << G4BestUnit((*(ener
463 << " e+ " << G4BestUnit((*(ener
464 << " proton " << G4BestUnit((*(ener
465 //// }
466 G4cout << G4endl;
467
468 if(aCouple->IsUsed())
469 {
470 G4cout << " Region(s) which use this cou
471 for(auto rItr=fG4RegionStore->cbegin();
472 rItr!=fG4RegionStore->cend(); +
473 {
474 if (IsCoupleUsedInTheRegion(aCouple, *
475 {
476 G4cout << " " << (*rItr)->GetName
477 }
478 }
479 }
480 }
481 G4cout << G4endl;
482 G4cout << "=================================
483 G4cout << G4endl;
484 }
485
486 // -------------------------------------------
487 G4bool G4ProductionCutsTable::StoreCutsTable(
488
489 {
490 // Store cuts and material information in fi
491
492 if (!StoreMaterialInfo(dir, ascii)) return f
493 if (!StoreMaterialCutsCoupleInfo(dir, ascii)
494 if (!StoreCutsInfo(dir, ascii)) return false
495
496 #ifdef G4VERBOSE
497 if (verboseLevel >2)
498 {
499 G4cout << "G4ProductionCutsTable::StoreCut
500 G4cout << " Material/Cuts information have
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::RetrieveCutsTab
517
518 {
519 if (!CheckForRetrieveCutsTable(dir, ascii))
520 if (!RetrieveCutsInfo(dir, ascii)) return fa
521 #ifdef G4VERBOSE
522 if (verboseLevel >2)
523 {
524 G4cout << "G4ProductionCutsTable::Retrieve
525 G4cout << " Material/Cuts information have
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::CheckForRetrieveCutsTab
543
544 {
545 // check stored material and cut values are
546 // with the current detector setup
547
548 G4cerr << "G4ProductionCutsTable::CheckForRe
549 // isNeedForRestoreCoupleInfo = false;
550 if (!CheckMaterialInfo(directory, ascii)) re
551 if (verboseLevel >2)
552 {
553 G4cerr << "G4ProductionCutsTable::CheckM
554 }
555 if (!CheckMaterialCutsCoupleInfo(directory,
556 if (verboseLevel >2)
557 {
558 G4cerr << "G4ProductionCutsTable::CheckMat
559 << G4endl;
560 }
561 return true;
562 }
563
564 // -------------------------------------------
565 G4bool G4ProductionCutsTable::StoreMaterialInf
566
567 {
568 // Store material information in files under
569
570 const G4String fileName = directory + "/" +
571 const G4String key = "MATERIAL-V3.0";
572 std::ofstream fOut;
573
574 // open output file
575 if (!ascii ) fOut.open(fileName,std::ios::o
576 else fOut.open(fileName,std::ios::o
577
578 // check if the file has been opened success
579 if (!fOut)
580 {
581 #ifdef G4VERBOSE
582 if (verboseLevel>0)
583 {
584 G4cerr << "G4ProductionCutsTable::StoreM
585 G4cerr << "Cannot open file: " << fileNa
586 }
587 #endif
588 G4Exception( "G4ProductionCutsTable::Store
589 "ProcCuts102", JustWarning, "
590 return false;
591 }
592
593 const G4MaterialTable* matTable = G4Material
594 // number of materials in the table
595 G4int numberOfMaterial = (G4int)matTable->si
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>
610 {
611 fOut << std::setw(FixedStringLengthForSt
612 << ((*matTable)[idx])->GetName();
613 fOut << std::setw(FixedStringLengthForSt
614 << ((*matTable)[idx])->GetDensity()
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<FixedStringL
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), si
639
640 // material name and density
641 for (std::size_t imat=0; static_cast<G4int
642 {
643 G4String name = ((*matTable)[imat])->Ge
644 G4double density = ((*matTable)[imat])->
645 for (i=0; i<FixedStringLengthForStore; +
646 temp[i] = '\0';
647 for (i=0; i<name.length() && i<FixedStri
648 temp[i]=name[(G4int)i];
649 fOut.write(temp, FixedStringLengthForSto
650 fOut.write( (char*)(&density), sizeof(G4
651 }
652 }
653
654 fOut.close();
655 return true;
656 }
657
658 // -------------------------------------------
659 G4bool G4ProductionCutsTable::CheckMaterialInf
660
661 {
662 // Check stored material is consistent with
663
664 const G4String fileName = directory + "/" +
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|
670 else fIn.open(fileName,std::ios::in)
671
672 // check if the file has been opened success
673 if (!fIn)
674 {
675 #ifdef G4VERBOSE
676 if (verboseLevel >0)
677 {
678 G4cerr << "G4ProductionCutsTable::CheckM
679 G4cerr << "Cannot open file: " << fileNa
680 }
681 #endif
682 G4Exception( "G4ProductionCutsTable::Check
683 "ProcCuts102", JustWarning, "
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::CheckM
706 G4cerr << "Key word in " << fileName <<
707 G4cerr <<"( should be "<< key << ")" <
708 }
709 #endif
710 G4Exception( "G4ProductionCutsTable::Check
711 "ProcCuts103", JustWarning, "
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::Check
728 "ProcCuts108", JustWarning,
729 "Number of materials is less
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::Chec
743 G4cout << "Encountered End of File " ;
744 G4cout << " at " << idx+1 << "th mate
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(G4dou
764 }
765 if (fIn.fail())
766 {
767 #ifdef G4VERBOSE
768 if (verboseLevel >0)
769 {
770 G4cerr << "G4ProductionCutsTable::Chec
771 G4cerr << "Bad data format ";
772 G4cerr << " at " << idx+1 << "th mate
773 }
774 #endif
775 G4Exception( "G4ProductionCutsTable::Che
776 "ProcCuts103", JustWarning,
777 fIn.close();
778 return false;
779 }
780
781 G4Material* aMaterial = G4Material::GetMat
782 if (aMaterial == nullptr ) continue;
783
784 G4double ratio = std::fabs(density/aMateri
785 if ((0.999>ratio) || (ratio>1.001) )
786 {
787 #ifdef G4VERBOSE
788 if (verboseLevel >0)
789 {
790 G4cerr << "G4ProductionCutsTable::Chec
791 G4cerr << "Inconsistent material densi
792 G4cerr << " at " << idx+1 << "th mate
793 G4cerr << "Name: " << name << G4endl
794 G4cerr << "Density:" << std::setiosfla
795 << density / (g/cm3) ;
796 G4cerr << "(should be " << aMaterial->
797 << " [g/cm3]"<< G4endl;
798 G4cerr << std::resetiosflags(std::ios:
799 }
800 #endif
801 G4Exception( "G4ProductionCutsTable::Che
802 "ProcCuts104", JustWarning,
803 fIn.close();
804 return false;
805 }
806 }
807
808 fIn.close();
809 return true;
810 }
811
812 // -------------------------------------------
813 G4bool
814 G4ProductionCutsTable::StoreMaterialCutsCouple
815
816 {
817 // Store materialCutsCouple information in f
818
819 const G4String fileName = directory + "/" +
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::ou
826 else fOut.open(fileName,std::ios::ou
827
828
829 // check if the file has been opened success
830 if (!fOut)
831 {
832 #ifdef G4VERBOSE
833 if (verboseLevel >0)
834 {
835 G4cerr << "G4ProductionCutsTable::StoreM
836 G4cerr << "Cannot open file: " << fileNa
837 }
838 #endif
839 G4Exception( "G4ProductionCutsTable::Store
840 "ProcCuts102",
841 JustWarning, "Cannot open fil
842 return false;
843 }
844 G4int numberOfCouples = (G4int)coupleTable.s
845 if (ascii)
846 {
847 /////////////// ASCII mode //////////////
848 // key word
849 fOut << std::setw(FixedStringLengthForStor
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<FixedStringL
862 temp[i]=key[(G4int)i];
863 fOut.write(temp, FixedStringLengthForStore
864
865 // number of couples in the table
866 fOut.write( (char*)(&numberOfCouples), siz
867 }
868
869 // Loop over all couples
870 for (auto cItr=coupleTable.cbegin(); cItr!=c
871 {
872 G4MaterialCutsCouple* aCouple = (*cItr);
873 G4int index = aCouple->GetIndex();
874 // cut value
875 G4ProductionCuts* aCut = aCouple->GetProdu
876 G4double cutValues[NumberOfG4CutIndex];
877 for (std::size_t idx=0; idx <NumberOfG4Cut
878 {
879 cutValues[idx] = aCut->GetProductionCut(
880 }
881 // material/region info
882 G4String materialName = aCouple->GetMateri
883 G4String regionName = "NONE";
884 if (aCouple->IsUsed())
885 {
886 for(auto rItr=fG4RegionStore->cbegin();
887 rItr!=fG4RegionStore->cend(); +
888 {
889 if (IsCoupleUsedInTheRegion(aCouple, *
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(FixedStringLengthForSt
905
906 // region name
907 fOut << std::setw(FixedStringLengthForSt
908
909 fOut.setf(std::ios::scientific);
910 // cut values
911 for (std::size_t idx=0; idx< NumberOfG4C
912 {
913 fOut << std::setw(FixedStringLengthFor
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(G4in
924
925 // material name
926 std::size_t i;
927 for (i=0; i<FixedStringLengthForStore; +
928 temp[i] = '\0';
929 for (i=0; i<materialName.length() && i<F
930 temp[i]=materialName[(G4int)i];
931 fOut.write(temp, FixedStringLengthForSto
932
933 // region name
934 for (i=0; i<FixedStringLengthForStore; +
935 temp[i] = '\0';
936 for (i=0; i<regionName.length() && i<Fix
937 temp[i]=regionName[(G4int)i];
938 fOut.write(temp, FixedStringLengthForSto
939
940 // cut values
941 for (std::size_t idx=0; idx< NumberOfG4C
942 {
943 fOut.write( (char*)(&(cutValues[idx])
944 }
945 }
946 }
947 fOut.close();
948 return true;
949 }
950
951 // -------------------------------------------
952 G4bool
953 G4ProductionCutsTable::CheckMaterialCutsCouple
954
955 {
956 // Check stored materialCutsCouple is consis
957 // with the current detector setup.
958
959 const G4String fileName = directory + "/" +
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|
965 else fIn.open(fileName,std::ios::in)
966
967 // check if the file has been opened success
968 if (!fIn)
969 {
970 #ifdef G4VERBOSE
971 if (verboseLevel >0)
972 {
973 G4cerr << "G4ProductionCutTable::CheckMa
974 G4cerr << "Cannot open file!" << fileNam
975 }
976 #endif
977 G4Exception( "G4ProductionCutsTable::Check
978 "ProcCuts102", JustWarning, "
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::CheckM
1001 G4cerr << "Key word in " << fileName <<
1002 G4cerr <<"( should be "<< key << ")"
1003 }
1004 #endif
1005 G4Exception( "G4ProductionCutsTable::Chec
1006 "ProcCuts103", JustWarning,
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), size
1020 }
1021
1022 // Reset MCCIndexConversionTable
1023 mccConversionTable.Reset(numberOfCouples);
1024
1025 // Read in couple information
1026 for (G4int idx=0; idx<numberOfCouples; ++id
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, FixedStringLengthFor
1047 }
1048 // read in index and region name
1049 char region_name[FixedStringLengthForStor
1050 if (ascii)
1051 {
1052 fIn >> region_name;
1053 }
1054 else
1055 {
1056 fIn.read(region_name, FixedStringLength
1057 }
1058 // cut value
1059 G4double cutValues[NumberOfG4CutIndex];
1060 for (std::size_t i=0; i< NumberOfG4CutInd
1061 {
1062 if (ascii)
1063 {
1064 fIn >> cutValues[i];
1065 cutValues[i] *= (mm);
1066 }
1067 else
1068 {
1069 fIn.read( (char*)(&(cutValues[i])), s
1070 }
1071 }
1072
1073 // Loop over all couples
1074 G4bool fOK = false;
1075 G4MaterialCutsCouple* aCouple = nullptr;
1076 for (auto cItr=coupleTable.cbegin(); cItr
1077 {
1078 aCouple = (*cItr);
1079 // check material name
1080 if ( mat_name != aCouple->GetMaterial(
1081 // check cut values
1082 G4ProductionCuts* aCut = aCouple->GetPr
1083 G4bool fRatio = true;
1084 for (std::size_t j=0; j< NumberOfG4CutI
1085 {
1086 // check ratio only if values are not
1087 if (cutValues[j] != aCut->GetProducti
1088 {
1089 G4double ratio = cutValues[j]/aCut
1090 fRatio = fRatio && (0.999<ratio) &&
1091 }
1092 }
1093 if (!fRatio) continue;
1094 // MCC matched
1095 fOK = true;
1096 mccConversionTable.SetNewIndex(index, a
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
1111 if (fRegion == nullptr)
1112 {
1113 G4cout << "G4ProductionCutTable::
1114 G4cout << "Region " << regionname
1115 G4cout << index << ": in " << fil
1116 }
1117 }
1118 if (((regionname == "NONE") && (aCou
1119 || ((fRegion!=nullptr) && !IsCoupleU
1120 {
1121 G4cout << "G4ProductionCutTable::Ch
1122 << G4endl;
1123 G4cout << "A Couple is used differe
1124 G4cout << index << ": in " << fileN
1125 G4cout << " material: " << mat_name
1126 G4cout << " region: " << region_nam
1127 for (std::size_t ii=0; ii< NumberOf
1128 {
1129 G4cout << "cut[" << ii << "]=" <<
1130 G4cout << " mm : ";
1131 }
1132 G4cout << G4endl;
1133 }
1134 else if ( index != aCouple->GetIndex
1135 {
1136 G4cout << "G4ProductionCutTable::Ch
1137 G4cout << "Index of couples was mod
1138 G4cout << aCouple->GetIndex() << ":
1139 << aCouple->GetMaterial()->G
1140 G4cout <<" is defined as " ;
1141 G4cout << index << ":" << mat_name
1142 }
1143 else
1144 {
1145 G4cout << "G4ProductionCutTable::Ch
1146 G4cout << index << ":" << mat_name
1147 G4cout << " is consistent with curr
1148 }
1149 }
1150 #endif
1151 }
1152 else
1153 {
1154 #ifdef G4VERBOSE
1155 if (verboseLevel >0)
1156 {
1157 G4cout << "G4ProductionCutTable::Chec
1158 << G4endl;
1159 G4cout << "Couples are not defined in
1160 G4cout << index << ": in " << fileNam
1161 G4cout << " material: " << mat_name ;
1162 G4cout << " region: " << region_name
1163 for (std::size_t ii=0; ii< NumberOfG4
1164 {
1165 G4cout << "cut[" << ii << "]=" << c
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(c
1179 G
1180 {
1181 // Store cut values information in files un
1182
1183 const G4String fileName = directory + "/" +
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::o
1190 else fOut.open(fileName,std::ios::o
1191
1192 // check if the file has been opened succes
1193 if (!fOut)
1194 {
1195 if(verboseLevel>0)
1196 {
1197 G4cerr << "G4ProductionCutsTable::Store
1198 G4cerr << "Cannot open file: " << fileN
1199 }
1200 G4Exception( "G4ProductionCutsTable::Stor
1201 "ProcCuts102", JustWarning,
1202 return false;
1203 }
1204
1205 G4int numberOfCouples = (G4int)coupleTable.
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; ++
1221 temp[i] = '\0';
1222 for (i=0; i<key.length() && i<FixedString
1223 temp[i]=key[(G4int)i];
1224 fOut.write(temp, FixedStringLengthForStor
1225
1226 // number of couples in the table
1227 fOut.write( (char*)(&numberOfCouples), si
1228 }
1229
1230 for (std::size_t idx=0; idx <NumberOfG4CutI
1231 {
1232 const std::vector<G4double>* fRange = Ge
1233 const std::vector<G4double>* fEnergy = Ge
1234 std::size_t i=0;
1235 // Loop over all couples
1236 for (auto cItr=coupleTable.cbegin();cItr!
1237 {
1238 if (ascii)
1239 {
1240 /////////////// ASCII mode /////////
1241 fOut.setf(std::ios::scientific);
1242 fOut << std::setw(20) << (*fRange)[i]
1243 fOut << std::setw(20) << (*fEnergy)[i
1244 fOut.unsetf(std::ios::scientific);
1245 }
1246 else
1247 {
1248 /////////////// Binary mode ////////
1249 G4double cut = (*fRange)[i];
1250 fOut.write((char*)(&cut), sizeof(G4do
1251 cut = (*fEnergy)[i];
1252 fOut.write((char*)(&cut), sizeof(G4do
1253 }
1254 }
1255 }
1256 fOut.close();
1257 return true;
1258 }
1259
1260 // ------------------------------------------
1261 G4bool G4ProductionCutsTable::RetrieveCutsInf
1262
1263 {
1264 // Retrieve cut values information in files
1265
1266 const G4String fileName = directory + "/" +
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
1272 else fIn.open(fileName,std::ios::in
1273
1274 // check if the file has been opened succes
1275 if (!fIn)
1276 {
1277 if (verboseLevel >0)
1278 {
1279 G4cerr << "G4ProductionCutTable::Retrie
1280 G4cerr << "Cannot open file: " << fileN
1281 }
1282 G4Exception( "G4ProductionCutsTable::Retr
1283 "ProcCuts102", JustWarning,
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::Retrie
1305 G4cerr << "Key word in " << fileName <<
1306 G4cerr <<"( should be "<< key << ")"
1307 }
1308 G4Exception( "G4ProductionCutsTable::Retr
1309 "ProcCuts103", JustWarning,
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::Re
1321 "ProcCuts103", JustWarning
1322 return false;
1323 }
1324 }
1325 else
1326 {
1327 fIn.read( (char*)(&numberOfCouples), size
1328 }
1329
1330 if (numberOfCouples > static_cast<G4int>(mc
1331 {
1332 G4Exception( "G4ProductionCutsTable::Retr
1333 "ProcCuts109", JustWarning,
1334 "Number of Couples in the fi
1335 }
1336 numberOfCouples = (G4int)mccConversionTable
1337
1338 for (std::size_t idx=0; static_cast<G4int>(
1339 {
1340 std::vector<G4double>* fRange = rangeCut
1341 std::vector<G4double>* fEnergy = energyCu
1342 fRange->clear();
1343 fEnergy->clear();
1344
1345 // Loop over all couples
1346 for (std::size_t i=0; static_cast<G4int>(
1347 {
1348 G4double rcut, ecut;
1349 if (ascii)
1350 {
1351 fIn >> rcut >> ecut;
1352 if (fIn.fail())
1353 {
1354 G4Exception( "G4ProductionCutsTable
1355 "ProcCuts103", JustWar
1356 return false;
1357 }
1358 rcut *= mm;
1359 ecut *= keV;
1360 }
1361 else
1362 {
1363 fIn.read((char*)(&rcut), sizeof(G4dou
1364 fIn.read((char*)(&ecut), sizeof(G4dou
1365 }
1366 if (!mccConversionTable.IsUsed(i)) cont
1367 std::size_t new_index = mccConversionTa
1368 (*fRange)[new_index] = rcut;
1369 (*fEnergy)[new_index] = ecut;
1370 }
1371 }
1372 return true;
1373 }
1374
1375 // ------------------------------------------
1376 void G4ProductionCutsTable::SetVerboseLevel(G
1377 {
1378 // Set same verbosity to all registered Ran
1379
1380 verboseLevel = value;
1381 for (G4int ip=0; ip< NumberOfG4CutIndex; ++
1382 {
1383 if (converters[ip] != nullptr )
1384 {
1385 converters[ip]->SetVerboseLevel(value);
1386 }
1387 }
1388 }
1389
1390 // ------------------------------------------
1391 G4double G4ProductionCutsTable::GetMaxEnergyC
1392 {
1393 return G4VRangeToEnergyConverter::GetMaxEne
1394 }
1395
1396 // ------------------------------------------
1397 void G4ProductionCutsTable::SetMaxEnergyCut(G
1398 {
1399 G4VRangeToEnergyConverter::SetMaxEnergyCut(
1400 }
1401