| 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
34 #include "BrachyUserScoreWriter.hh" 34 #include "BrachyUserScoreWriter.hh"
35 #include "G4AnalysisManager.hh" 35 #include "G4AnalysisManager.hh"
36 #include "G4MultiFunctionalDetector.hh" 36 #include "G4MultiFunctionalDetector.hh"
37 #include "G4SDParticleFilter.hh" 37 #include "G4SDParticleFilter.hh"
38 #include "G4VPrimitiveScorer.hh" 38 #include "G4VPrimitiveScorer.hh"
39 #include "G4VScoringMesh.hh" 39 #include "G4VScoringMesh.hh"
40 #include "G4SystemOfUnits.hh" 40 #include "G4SystemOfUnits.hh"
41 #include <map> 41 #include <map>
42 #include <fstream> 42 #include <fstream>
43 // The default output is 43 // The default output is
44 // voxelX, voxelY, voxelZ, edep 44 // voxelX, voxelY, voxelZ, edep
45 // The BrachyUserScoreWriter allows to change 45 // The BrachyUserScoreWriter allows to change the format of the output file.
46 // in the specific case: 46 // in the specific case:
47 // xx (mm) yy(mm) zz(mm) edep(keV) 47 // xx (mm) yy(mm) zz(mm) edep(keV)
48 // The same information is stored in a ntuple, 48 // The same information is stored in a ntuple, in the
49 // brachytherapy.root file 49 // brachytherapy.root file
50 50
51 BrachyUserScoreWriter::BrachyUserScoreWriter() 51 BrachyUserScoreWriter::BrachyUserScoreWriter():
52 G4VScoreWriter() 52 G4VScoreWriter()
53 { 53 {
54 } 54 }
55 55
56 BrachyUserScoreWriter::~BrachyUserScoreWriter( 56 BrachyUserScoreWriter::~BrachyUserScoreWriter()
57 {;} 57 {;}
58 58
59 void BrachyUserScoreWriter::DumpQuantityToFile 59 void BrachyUserScoreWriter::DumpQuantityToFile(const G4String & psName,
60 60 const G4String & fileName,
61 61 const G4String & option)
62 { 62 {
63 using MeshScoreMap = G4VScoringMesh::MeshScore 63 using MeshScoreMap = G4VScoringMesh::MeshScoreMap;
64 64
65 if(verboseLevel > 0) 65 if(verboseLevel > 0)
66 {G4cout << "BrachyUserScorer-defined DumpQua 66 {G4cout << "BrachyUserScorer-defined DumpQuantityToFile() method is invoked."
67 << G4endl; 67 << G4endl;
68 } 68 }
69 69
70 // change the option string into lowercase to 70 // change the option string into lowercase to the case-insensitive.
71 G4String opt = option; 71 G4String opt = option;
72 std::transform(opt.begin(), opt.end(), opt.beg 72 std::transform(opt.begin(), opt.end(), opt.begin(), (int (*)(int))(tolower));
73 73
74 // confirm the option 74 // confirm the option
75 if(opt.size() == 0) opt = "csv"; 75 if(opt.size() == 0) opt = "csv";
76 76
77 // open the file 77 // open the file
78 std::ofstream ofile(fileName); 78 std::ofstream ofile(fileName);
79 79
80 if(!ofile) 80 if(!ofile)
81 { 81 {
82 G4cerr << "ERROR : DumpToFile : File open e 82 G4cerr << "ERROR : DumpToFile : File open error -> " << fileName << G4endl;
83 return; 83 return;
84 } 84 }
85 ofile << "# mesh name: " << fScoringMesh->Ge 85 ofile << "# mesh name: " << fScoringMesh->GetWorldName() << G4endl;
86 86
87 // retrieve the map 87 // retrieve the map
88 MeshScoreMap fSMap = fScoringMesh -> GetScoreM 88 MeshScoreMap fSMap = fScoringMesh -> GetScoreMap();
89 89
90 auto msMapItr = fSMap.find(psName); 90 auto msMapItr = fSMap.find(psName);
91 91
92 if(msMapItr == fSMap.end()) 92 if(msMapItr == fSMap.end())
93 { 93 {
94 G4cerr << "ERROR : DumpToFile : Unknown qua 94 G4cerr << "ERROR : DumpToFile : Unknown quantity, \""<< psName
95 << "\"." << G4endl; 95 << "\"." << G4endl;
96 return; 96 return;
97 } 97 }
98 98
99 auto score = msMapItr-> second-> GetMap(); 99 auto score = msMapItr-> second-> GetMap();
100 100
101 ofile << "# primitive scorer name: " << msMapI 101 ofile << "# primitive scorer name: " << msMapItr -> first << G4endl;
102 // 102 //
103 // Write quantity in the ASCII output file and 103 // Write quantity in the ASCII output file and in brachytherapy.root
104 // 104 //
105 ofile << std::setprecision(16); // for double 105 ofile << std::setprecision(16); // for double value with 8 bytes
106 106
107 auto analysisManager = G4AnalysisManager::Inst 107 auto analysisManager = G4AnalysisManager::Instance();
108 108
109 G4bool fileOpen = analysisManager -> OpenFile( 109 G4bool fileOpen = analysisManager -> OpenFile("brachytherapy.root");
110 if (! fileOpen) { 110 if (! fileOpen) {
111 G4cerr << "\n---> The ROOT output file has 111 G4cerr << "\n---> The ROOT output file has not been opened "
112 << analysisManager->GetFileName() < 112 << analysisManager->GetFileName() << G4endl;
113 } 113 }
114 114
115 G4cout << "Using " << analysisManager -> GetTy 115 G4cout << "Using " << analysisManager -> GetType() << G4endl;
116 analysisManager -> SetVerboseLevel(1); 116 analysisManager -> SetVerboseLevel(1);
117 analysisManager -> SetActivation(true); 117 analysisManager -> SetActivation(true);
118 118
119 // Create histograms 119 // Create histograms
120 G4int histo2= analysisManager-> CreateH2("h20" 120 G4int histo2= analysisManager-> CreateH2("h20","edep2Dxy", 801, -100.125, 100.125, 801, -100.125, 100.125);
121 121
122 // Histo 0 with the energy spectrum will not b 122 // Histo 0 with the energy spectrum will not be saved
123 // in brachytherapy.root 123 // in brachytherapy.root
124 analysisManager->SetH1Activation(0, false); 124 analysisManager->SetH1Activation(0, false);
125 analysisManager->SetH2Activation(histo2, true) 125 analysisManager->SetH2Activation(histo2, true);
126 126
127 for(int x = 0; x < fNMeshSegments[0]; x++) { 127 for(int x = 0; x < fNMeshSegments[0]; x++) {
128 for(int y = 0; y < fNMeshSegments[1]; y++) 128 for(int y = 0; y < fNMeshSegments[1]; y++) {
129 for(int z = 0; z < fNMeshSegments[2]; z++ 129 for(int z = 0; z < fNMeshSegments[2]; z++){
130 G4int numberOfVoxel_x = fNMeshSegments 130 G4int numberOfVoxel_x = fNMeshSegments[0];
131 G4int numberOfVoxel_y = fNMeshSegments 131 G4int numberOfVoxel_y = fNMeshSegments[1];
132 G4int numberOfVoxel_z =fNMeshSegments[ 132 G4int numberOfVoxel_z =fNMeshSegments[2];
133 // If the voxel width is changed in th 133 // If the voxel width is changed in the macro file,
134 // the voxel width variable must be up 134 // the voxel width variable must be updated
135 G4double voxelWidth = 0.25 *CLHEP::mm; 135 G4double voxelWidth = 0.25 *CLHEP::mm;
136 // 136 //
137 G4double xx = ( - numberOfVoxel_x + 1+ 137 G4double xx = ( - numberOfVoxel_x + 1+ 2*x )* voxelWidth/2;
138 G4double yy = ( - numberOfVoxel_y + 1+ 138 G4double yy = ( - numberOfVoxel_y + 1+ 2*y )* voxelWidth/2;
139 G4double zz = ( - numberOfVoxel_z + 1+ 139 G4double zz = ( - numberOfVoxel_z + 1+ 2*z )* voxelWidth/2;
140 G4int idx = GetIndex(x, y, z); 140 G4int idx = GetIndex(x, y, z);
141 std::map<G4int, G4StatDouble*>::iterat 141 std::map<G4int, G4StatDouble*>::iterator value = score -> find(idx);
142 142
143 if (value != score -> end()) 143 if (value != score -> end())
144 { 144 {
145 // Print in the ASCII output file the 145 // Print in the ASCII output file the information
146 146
147 ofile << xx << " " << yy << " " << 147 ofile << xx << " " << yy << " " << zz <<" "
148 <<(value->second->sum_wx())/keV 148 <<(value->second->sum_wx())/keV << G4endl;
149 149
150 // Save the same information in the RO 150 // Save the same information in the ROOT output file
151 151
152 if(zz> -0.125 *CLHEP::mm && zz < 0.125/mm) 152 if(zz> -0.125 *CLHEP::mm && zz < 0.125/mm)
153 analysisManager->FillH2(histo2, xx, y 153 analysisManager->FillH2(histo2, xx, yy, (value->second->sum_wx())/keV);
154 }}}} 154 }}}}
155 155
156 ofile << std::setprecision(6); 156 ofile << std::setprecision(6);
157 157
158 // Close the output ASCII file 158 // Close the output ASCII file
159 ofile.close(); 159 ofile.close();
160 160
161 // Close the output ROOT file 161 // Close the output ROOT file
162 analysisManager -> Write(); 162 analysisManager -> Write();
163 analysisManager -> CloseFile(); 163 analysisManager -> CloseFile();
164 } 164 }
165 165