| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 /* 27 /*
28 // Code developed by: 28 // Code developed by:
29 // S.Guatelli, susanna@uow.edu.au 29 // S.Guatelli, susanna@uow.edu.au
30 // 30 //
31 Original code from geant4/examples/extended/ru 31 Original code from geant4/examples/extended/runAndEvent/RE03, by M. Asai
32 */ 32 */
33 << 33 #include <map>
>> 34 #include <fstream>
34 #include "BrachyUserScoreWriter.hh" 35 #include "BrachyUserScoreWriter.hh"
35 #include "G4AnalysisManager.hh" << 36 #include "G4SystemOfUnits.hh"
>> 37
>> 38 #ifdef ANALYSIS_USE
>> 39 #include "BrachyAnalysisManager.hh"
>> 40 #endif
>> 41
36 #include "G4MultiFunctionalDetector.hh" 42 #include "G4MultiFunctionalDetector.hh"
37 #include "G4SDParticleFilter.hh" 43 #include "G4SDParticleFilter.hh"
38 #include "G4VPrimitiveScorer.hh" 44 #include "G4VPrimitiveScorer.hh"
39 #include "G4VScoringMesh.hh" 45 #include "G4VScoringMesh.hh"
40 #include "G4SystemOfUnits.hh" <<
41 #include <map> <<
42 #include <fstream> <<
43 // The default output is 46 // The default output is
44 // voxelX, voxelY, voxelZ, edep 47 // voxelX, voxelY, voxelZ, edep
45 // The BrachyUserScoreWriter allows to change 48 // The BrachyUserScoreWriter allows to change the format of the output file.
46 // in the specific case: 49 // in the specific case:
47 // xx (mm) yy(mm) zz(mm) edep(keV) 50 // xx (mm) yy(mm) zz(mm) edep(keV)
48 // The same information is stored in a ntuple, 51 // The same information is stored in a ntuple, in the
49 // brachytherapy.root file 52 // brachytherapy.root file
50 53
51 BrachyUserScoreWriter::BrachyUserScoreWriter() 54 BrachyUserScoreWriter::BrachyUserScoreWriter():
52 G4VScoreWriter() 55 G4VScoreWriter()
53 { 56 {
54 } 57 }
55 58
56 BrachyUserScoreWriter::~BrachyUserScoreWriter( 59 BrachyUserScoreWriter::~BrachyUserScoreWriter()
57 {;} 60 {;}
58 61
59 void BrachyUserScoreWriter::DumpQuantityToFile 62 void BrachyUserScoreWriter::DumpQuantityToFile(const G4String & psName,
60 63 const G4String & fileName,
61 64 const G4String & option)
62 { 65 {
63 using MeshScoreMap = G4VScoringMesh::MeshScore <<
64 <<
65 if(verboseLevel > 0) 66 if(verboseLevel > 0)
66 {G4cout << "BrachyUserScorer-defined DumpQua 67 {G4cout << "BrachyUserScorer-defined DumpQuantityToFile() method is invoked."
67 << G4endl; 68 << G4endl;
68 } 69 }
69 70
70 // change the option string into lowercase to 71 // change the option string into lowercase to the case-insensitive.
71 G4String opt = option; 72 G4String opt = option;
72 std::transform(opt.begin(), opt.end(), opt.beg 73 std::transform(opt.begin(), opt.end(), opt.begin(), (int (*)(int))(tolower));
73 74
74 // confirm the option 75 // confirm the option
75 if(opt.size() == 0) opt = "csv"; 76 if(opt.size() == 0) opt = "csv";
76 77
77 // open the file 78 // open the file
78 std::ofstream ofile(fileName); 79 std::ofstream ofile(fileName);
79 80
80 if(!ofile) 81 if(!ofile)
81 { 82 {
82 G4cerr << "ERROR : DumpToFile : File open e 83 G4cerr << "ERROR : DumpToFile : File open error -> " << fileName << G4endl;
83 return; 84 return;
84 } 85 }
85 ofile << "# mesh name: " << fScoringMesh->Ge 86 ofile << "# mesh name: " << fScoringMesh->GetWorldName() << G4endl;
86 87
87 // retrieve the map 88 // retrieve the map
88 MeshScoreMap fSMap = fScoringMesh -> GetScoreM 89 MeshScoreMap fSMap = fScoringMesh -> GetScoreMap();
89 << 90
90 auto msMapItr = fSMap.find(psName); << 91 MeshScoreMap::const_iterator msMapItr = fSMap.find(psName);
91 92
92 if(msMapItr == fSMap.end()) 93 if(msMapItr == fSMap.end())
93 { 94 {
94 G4cerr << "ERROR : DumpToFile : Unknown qua 95 G4cerr << "ERROR : DumpToFile : Unknown quantity, \""<< psName
95 << "\"." << G4endl; 96 << "\"." << G4endl;
96 return; 97 return;
97 } 98 }
98 99
99 auto score = msMapItr-> second-> GetMap(); << 100 std::map<G4int, G4StatDouble*> * score = msMapItr -> second-> GetMap();
100 101
101 ofile << "# primitive scorer name: " << msMapI 102 ofile << "# primitive scorer name: " << msMapItr -> first << G4endl;
102 // 103 //
103 // Write quantity in the ASCII output file and 104 // Write quantity in the ASCII output file and in brachytherapy.root
104 // 105 //
105 ofile << std::setprecision(16); // for double 106 ofile << std::setprecision(16); // for double value with 8 bytes
106 <<
107 auto analysisManager = G4AnalysisManager::Inst <<
108 <<
109 G4bool fileOpen = analysisManager -> OpenFile( <<
110 if (! fileOpen) { <<
111 G4cerr << "\n---> The ROOT output file has <<
112 << analysisManager->GetFileName() < <<
113 } <<
114 <<
115 G4cout << "Using " << analysisManager -> GetTy <<
116 analysisManager -> SetVerboseLevel(1); <<
117 analysisManager -> SetActivation(true); <<
118 <<
119 // Create histograms <<
120 G4int histo2= analysisManager-> CreateH2("h20" <<
121 <<
122 // Histo 0 with the energy spectrum will not b <<
123 // in brachytherapy.root <<
124 analysisManager->SetH1Activation(0, false); <<
125 analysisManager->SetH2Activation(histo2, true) <<
126 107
127 for(int x = 0; x < fNMeshSegments[0]; x++) { 108 for(int x = 0; x < fNMeshSegments[0]; x++) {
128 for(int y = 0; y < fNMeshSegments[1]; y++) 109 for(int y = 0; y < fNMeshSegments[1]; y++) {
129 for(int z = 0; z < fNMeshSegments[2]; z++ 110 for(int z = 0; z < fNMeshSegments[2]; z++){
130 G4int numberOfVoxel_x = fNMeshSegments 111 G4int numberOfVoxel_x = fNMeshSegments[0];
131 G4int numberOfVoxel_y = fNMeshSegments 112 G4int numberOfVoxel_y = fNMeshSegments[1];
132 G4int numberOfVoxel_z =fNMeshSegments[ 113 G4int numberOfVoxel_z =fNMeshSegments[2];
133 // If the voxel width is changed in th 114 // If the voxel width is changed in the macro file,
134 // the voxel width variable must be up 115 // the voxel width variable must be updated
135 G4double voxelWidth = 0.25 *CLHEP::mm; << 116 G4double voxelWidth = 0.25 *mm;
136 // 117 //
137 G4double xx = ( - numberOfVoxel_x + 1+ 118 G4double xx = ( - numberOfVoxel_x + 1+ 2*x )* voxelWidth/2;
138 G4double yy = ( - numberOfVoxel_y + 1+ 119 G4double yy = ( - numberOfVoxel_y + 1+ 2*y )* voxelWidth/2;
139 G4double zz = ( - numberOfVoxel_z + 1+ 120 G4double zz = ( - numberOfVoxel_z + 1+ 2*z )* voxelWidth/2;
140 G4int idx = GetIndex(x, y, z); 121 G4int idx = GetIndex(x, y, z);
141 std::map<G4int, G4StatDouble*>::iterat 122 std::map<G4int, G4StatDouble*>::iterator value = score -> find(idx);
142 << 123
143 if (value != score -> end()) 124 if (value != score -> end())
144 { 125 {
145 // Print in the ASCII output file the 126 // Print in the ASCII output file the information
146 127
147 ofile << xx << " " << yy << " " << 128 ofile << xx << " " << yy << " " << zz <<" "
148 <<(value->second->sum_wx())/keV 129 <<(value->second->sum_wx())/keV << G4endl;
149 << 130
150 // Save the same information in the RO << 131 #ifdef ANALYSIS_USE
>> 132 // Save the same information in the output analysis file
>> 133 BrachyAnalysisManager* analysis = BrachyAnalysisManager::GetInstance();
151 134
152 if(zz> -0.125 *CLHEP::mm && zz < 0.125/mm) << 135 if(zz> -0.125 *mm && zz < 0.125*mm) analysis -> FillH2WithEnergyDeposition(xx,yy, (value->second->sum_wx())/keV);
153 analysisManager->FillH2(histo2, xx, y << 136 #endif
154 }}}} 137 }}}}
155 138
156 ofile << std::setprecision(6); 139 ofile << std::setprecision(6);
157 140
158 // Close the output ASCII file 141 // Close the output ASCII file
159 ofile.close(); 142 ofile.close();
160 <<
161 // Close the output ROOT file <<
162 analysisManager -> Write(); <<
163 analysisManager -> CloseFile(); <<
164 } 143 }
165 144