| 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 // $Id: XrayTelAnalysis.cc 68710 2013-04-05 09:04:21Z gcosmo $
27 // 28 //
28 // Author: A. Pfeiffer (Andreas.Pfeiffer@cern 29 // Author: A. Pfeiffer (Andreas.Pfeiffer@cern.ch)
29 // (copied from his UserAnalyser class 30 // (copied from his UserAnalyser class)
30 // 31 //
31 // History: 32 // History:
32 // ----------- 33 // -----------
33 // 19 Mar 2013 LP Migrated to G4tool 34 // 19 Mar 2013 LP Migrated to G4tools
34 // 8 Nov 2002 GS Migration to AIDA 35 // 8 Nov 2002 GS Migration to AIDA 3
35 // 7 Nov 2001 MGP Implementation 36 // 7 Nov 2001 MGP Implementation
36 // 37 //
37 // ------------------------------------------- 38 // -------------------------------------------------------------------
38 39
39 #include <fstream> 40 #include <fstream>
40 #include <iomanip> 41 #include <iomanip>
41 #include "G4AutoLock.hh" << 42
42 #include "XrayTelAnalysis.hh" 43 #include "XrayTelAnalysis.hh"
43 #include "globals.hh" 44 #include "globals.hh"
44 #include "G4SystemOfUnits.hh" 45 #include "G4SystemOfUnits.hh"
45 #include "G4Track.hh" 46 #include "G4Track.hh"
46 #include "G4ios.hh" 47 #include "G4ios.hh"
47 #include "G4SteppingManager.hh" 48 #include "G4SteppingManager.hh"
48 #include "G4ThreeVector.hh" 49 #include "G4ThreeVector.hh"
49 50
50 XrayTelAnalysis* XrayTelAnalysis::instance = 0 51 XrayTelAnalysis* XrayTelAnalysis::instance = 0;
51 52
52 namespace { <<
53 //Mutex to acquire access to singleton insta <<
54 G4Mutex instanceMutex = G4MUTEX_INITIALIZER; <<
55 //Mutex to acquire accss to histograms creat <<
56 //It is also used to control all operations <<
57 //File writing and check analysis <<
58 G4Mutex dataManipulationMutex = G4MUTEX_INIT <<
59 } <<
60 <<
61 XrayTelAnalysis::XrayTelAnalysis() : 53 XrayTelAnalysis::XrayTelAnalysis() :
62 asciiFile(0),nEnteringTracks(0), totEntering << 54 asciiFile(0)
63 { 55 {
64 histFileName = "xraytel"; 56 histFileName = "xraytel";
65 57
66 58
67 asciiFileName="xraytel.out"; 59 asciiFileName="xraytel.out";
68 asciiFile = new std::ofstream(asciiFileName) 60 asciiFile = new std::ofstream(asciiFileName);
69 61
70 if(asciiFile->is_open()) 62 if(asciiFile->is_open())
71 (*asciiFile) << "Energy (keV) x (mm) y 63 (*asciiFile) << "Energy (keV) x (mm) y (mm) z (mm)" << G4endl << G4endl;
72 } 64 }
73 65
74 XrayTelAnalysis::~XrayTelAnalysis() 66 XrayTelAnalysis::~XrayTelAnalysis()
75 { 67 {
76 if (asciiFile) 68 if (asciiFile)
77 delete asciiFile; 69 delete asciiFile;
78 if (nEnteringTracks) <<
79 delete nEnteringTracks; <<
80 if (totEnteringEnergy) <<
81 delete totEnteringEnergy; <<
82 } 70 }
83 71
84 72
85 XrayTelAnalysis* XrayTelAnalysis::getInstance( 73 XrayTelAnalysis* XrayTelAnalysis::getInstance()
86 { 74 {
87 G4AutoLock l(&instanceMutex); << 75 if (instance == 0) instance = new XrayTelAnalysis;
88 if (instance == 0) <<
89 instance = new XrayTelAnalysis; <<
90 return instance; 76 return instance;
91 } 77 }
92 78
93 79
94 void XrayTelAnalysis::book(G4bool isMaster) << 80 void XrayTelAnalysis::book()
95 { 81 {
96 G4AutoLock l(&dataManipulationMutex); <<
97 <<
98 //reset counters: do be done only once, by t <<
99 if (isMaster) <<
100 { <<
101 if (nEnteringTracks) <<
102 { <<
103 delete nEnteringTracks; <<
104 nEnteringTracks = 0; <<
105 } <<
106 nEnteringTracks = new std::map<G4int,G4i <<
107 <<
108 if (totEnteringEnergy) <<
109 { <<
110 delete totEnteringEnergy; <<
111 totEnteringEnergy = 0; <<
112 } <<
113 totEnteringEnergy = new std::map<G4int,G <<
114 } <<
115 <<
116 // Get/create analysis manager 82 // Get/create analysis manager
117 G4AnalysisManager* man = G4AnalysisManager:: 83 G4AnalysisManager* man = G4AnalysisManager::Instance();
118 man->SetDefaultFileType("root"); <<
119 84
120 // Open an output file: it is done in master << 85 // Open an output file
121 // printout is done only by the master, for << 86 G4cout << "Opening output file " << histFileName << " ... ";
122 if (isMaster) <<
123 G4cout << "Opening output file " << histFi <<
124 man->OpenFile(histFileName); 87 man->OpenFile(histFileName);
125 man->SetFirstHistoId(1); 88 man->SetFirstHistoId(1);
126 if (isMaster) << 89 G4cout << " done" << G4endl;
127 G4cout << " done" << G4endl; <<
128 90
129 // Book 1D histograms 91 // Book 1D histograms
130 man->CreateH1("h1","Energy, all /keV", 100, << 92 man->CreateH1("1","Energy, all /keV", 100,0.,100.);
131 man->CreateH1("h2","Energy, entering detecto << 93 man->CreateH1("2","Energy, entering detector /keV", 500,0.,500.);
132 94
133 // Book 2D histograms (notice: the numbering 95 // Book 2D histograms (notice: the numbering is independent)
134 man->CreateH2("d1","y-z, all /mm", 100,-500. << 96 man->CreateH2("1","y-z, all /mm", 100,-500.,500.,100,-500.,500.);
135 man->CreateH2("d2","y-z, entering detector / << 97 man->CreateH2("2","y-z, entering detector /mm", 200,-50.,50.,200,-50.,50.);
136 98
137 // Book ntuples 99 // Book ntuples
138 man->CreateNtuple("tree", "Track ntuple"); << 100 man->CreateNtuple("10", "Track ntuple");
139 man->CreateNtupleDColumn("energy"); 101 man->CreateNtupleDColumn("energy");
140 man->CreateNtupleDColumn("x"); 102 man->CreateNtupleDColumn("x");
141 man->CreateNtupleDColumn("y"); 103 man->CreateNtupleDColumn("y");
142 man->CreateNtupleDColumn("z"); 104 man->CreateNtupleDColumn("z");
143 man->CreateNtupleDColumn("dirx"); 105 man->CreateNtupleDColumn("dirx");
144 man->CreateNtupleDColumn("diry"); 106 man->CreateNtupleDColumn("diry");
145 man->CreateNtupleDColumn("dirz"); 107 man->CreateNtupleDColumn("dirz");
146 man->FinishNtuple(); 108 man->FinishNtuple();
147 } 109 }
148 110
149 111
150 void XrayTelAnalysis::finish(G4bool isMaster) << 112 void XrayTelAnalysis::finish()
151 { 113 {
152 G4AutoLock l(&dataManipulationMutex); << 114 asciiFile->close();
>> 115
153 // Save histograms 116 // Save histograms
154 G4AnalysisManager* man = G4AnalysisManager:: 117 G4AnalysisManager* man = G4AnalysisManager::Instance();
155 man->Write(); 118 man->Write();
156 man->CloseFile(); 119 man->CloseFile();
157 man->Clear(); << 120 // Complete clean-up
158 << 121 delete G4AnalysisManager::Instance();
159 if (!isMaster) <<
160 return; <<
161 <<
162 //only master performs these operations <<
163 asciiFile->close(); <<
164 <<
165 //Sequential run: output results! <<
166 if (nEnteringTracks->count(-2)) <<
167 { <<
168 G4cout << "End of Run summary (sequentia <<
169 G4cout << "Total Entering Detector : " < <<
170 << G4endl; <<
171 G4cout << "Total Entering Detector Energ <<
172 << (totEnteringEnergy->find(-2)->second <<
173 << " MeV" <<
174 << G4endl; <<
175 return; <<
176 } <<
177 <<
178 <<
179 //MT run: sum results <<
180 <<
181 <<
182 //MT build, but sequential run <<
183 if (nEnteringTracks->count(-1)) <<
184 { <<
185 G4cout << "End of Run summary (sequentia <<
186 G4cout << "Total Entering Detector : " < <<
187 << G4endl; <<
188 G4cout << "Total Entering Detector Energ <<
189 << (totEnteringEnergy->find(-1)->second <<
190 << " MeV" <<
191 << G4endl; <<
192 G4cout << "############################# <<
193 return; <<
194 } <<
195 <<
196 G4bool loopAgain = true; <<
197 G4int totEntries = 0; <<
198 G4double totEnergy = 0.; <<
199 <<
200 G4cout << "################################# <<
201 for (size_t i=0; loopAgain; i++) <<
202 { <<
203 //ok, this thread was found <<
204 if (nEnteringTracks->count(i)) <<
205 { <<
206 G4cout << "End of Run summary (thread= " < <<
207 G4int part = nEnteringTracks->find(i)->sec <<
208 G4cout << "Total Entering Detector : " << <<
209 G4double ene = totEnteringEnergy->find(i)- <<
210 G4cout << "Total Entering Detector Energy <<
211 << ene/MeV <<
212 << " MeV" << G4endl; <<
213 totEntries += part; <<
214 totEnergy += ene; <<
215 G4cout << "############################### <<
216 } <<
217 else <<
218 loopAgain = false; <<
219 } <<
220 //Report global numbers <<
221 if (totEntries) <<
222 { <<
223 G4cout << "End of Run summary (MT) globa <<
224 G4cout << "Total Entering Detector : " < <<
225 G4cout << "Total Entering Detector Energ <<
226 << totEnergy/MeV <<
227 << " MeV" << G4endl; <<
228 G4cout << "############################# <<
229 } <<
230 } 122 }
231 123
232 void XrayTelAnalysis::analyseStepping(const G4 124 void XrayTelAnalysis::analyseStepping(const G4Track& track, G4bool entering)
233 { 125 {
234 G4AutoLock l(&dataManipulationMutex); <<
235 eKin = track.GetKineticEnergy()/keV; 126 eKin = track.GetKineticEnergy()/keV;
236 G4ThreeVector pos = track.GetPosition()/mm; 127 G4ThreeVector pos = track.GetPosition()/mm;
237 y = pos.y(); 128 y = pos.y();
238 z = pos.z(); 129 z = pos.z();
239 G4ThreeVector dir = track.GetMomentumDirecti 130 G4ThreeVector dir = track.GetMomentumDirection();
240 dirX = dir.x(); 131 dirX = dir.x();
241 dirY = dir.y(); 132 dirY = dir.y();
242 dirZ = dir.z(); 133 dirZ = dir.z();
243 134
244 // Fill histograms 135 // Fill histograms
245 G4AnalysisManager* man = G4AnalysisManager:: 136 G4AnalysisManager* man = G4AnalysisManager::Instance();
246 man->FillH1(1,eKin); 137 man->FillH1(1,eKin);
247 man->FillH2(1,y,z); 138 man->FillH2(1,y,z);
248 139
249 // Fill histograms and ntuple, tracks enteri 140 // Fill histograms and ntuple, tracks entering the detector
250 if (entering) { 141 if (entering) {
251 // Fill and plot histograms 142 // Fill and plot histograms
252 man->FillH1(2,eKin); 143 man->FillH1(2,eKin);
253 man->FillH2(2,y,z); 144 man->FillH2(2,y,z);
254 145
255 man->FillNtupleDColumn(0,eKin); 146 man->FillNtupleDColumn(0,eKin);
256 man->FillNtupleDColumn(1,x); 147 man->FillNtupleDColumn(1,x);
257 man->FillNtupleDColumn(2,y); 148 man->FillNtupleDColumn(2,y);
258 man->FillNtupleDColumn(3,z); 149 man->FillNtupleDColumn(3,z);
259 man->FillNtupleDColumn(4,dirX); 150 man->FillNtupleDColumn(4,dirX);
260 man->FillNtupleDColumn(5,dirY); 151 man->FillNtupleDColumn(5,dirY);
261 man->FillNtupleDColumn(6,dirZ); 152 man->FillNtupleDColumn(6,dirZ);
262 man->AddNtupleRow(); 153 man->AddNtupleRow();
263 } 154 }
264 155
265 156
266 // Write to file 157 // Write to file
267 if (entering) { 158 if (entering) {
268 if(asciiFile->is_open()) { 159 if(asciiFile->is_open()) {
269 (*asciiFile) << std::setiosflags(std::io 160 (*asciiFile) << std::setiosflags(std::ios::fixed)
270 << std::setprecision(3) 161 << std::setprecision(3)
271 << std::setiosflags(std::ios::right) 162 << std::setiosflags(std::ios::right)
272 << std::setw(10); 163 << std::setw(10);
273 (*asciiFile) << eKin; 164 (*asciiFile) << eKin;
274 (*asciiFile) << std::setiosflags(std::io 165 (*asciiFile) << std::setiosflags(std::ios::fixed)
275 << std::setprecision(3) 166 << std::setprecision(3)
276 << std::setiosflags(std::ios::right) 167 << std::setiosflags(std::ios::right)
277 << std::setw(10); 168 << std::setw(10);
278 (*asciiFile) << x; 169 (*asciiFile) << x;
279 (*asciiFile) << std::setiosflags(std::io 170 (*asciiFile) << std::setiosflags(std::ios::fixed)
280 << std::setprecision(3) 171 << std::setprecision(3)
281 << std::setiosflags(std::ios::right) 172 << std::setiosflags(std::ios::right)
282 << std::setw(10); 173 << std::setw(10);
283 (*asciiFile) << y; 174 (*asciiFile) << y;
284 (*asciiFile) << std::setiosflags(std::io 175 (*asciiFile) << std::setiosflags(std::ios::fixed)
285 << std::setprecision(3) 176 << std::setprecision(3)
286 << std::setiosflags(std::ios::right) 177 << std::setiosflags(std::ios::right)
287 << std::setw(10); 178 << std::setw(10);
288 (*asciiFile) << z 179 (*asciiFile) << z
289 << G4endl; 180 << G4endl;
290 } 181 }
291 } 182 }
292 } <<
293 <<
294 void XrayTelAnalysis::Update(G4double energy,G <<
295 { <<
296 G4AutoLock l(&dataManipulationMutex); <<
297 //It already exists: increase the counter <<
298 if (nEnteringTracks->count(threadID)) <<
299 { <<
300 (nEnteringTracks->find(threadID)->second <<
301 } <<
302 else //enter a new one <<
303 { <<
304 G4int tracks = 1; <<
305 nEnteringTracks->insert(std::make_pair(t <<
306 } <<
307 <<
308 //It already exists: increase the counter <<
309 if (totEnteringEnergy->count(threadID)) <<
310 (totEnteringEnergy->find(threadID)->second <<
311 else //enter a new one <<
312 totEnteringEnergy->insert(std::make_pair(t <<
313 <<
314 return; <<
315 } 183 }
316 184
317 185