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