| 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 //
27 //
28 // Author: A. Pfeiffer (Andreas.Pfeiffer@cern
29 // (copied from his UserAnalyser class
30 //
31 // History:
32 // -----------
33 // 19 Mar 2013 LP Migrated to G4tool
34 // 8 Nov 2002 GS Migration to AIDA
35 // 7 Nov 2001 MGP Implementation
36 //
37 // -------------------------------------------
38
39 #include <fstream>
40 #include <iomanip>
41 #include "G4AutoLock.hh"
42 #include "XrayTelAnalysis.hh"
43 #include "globals.hh"
44 #include "G4SystemOfUnits.hh"
45 #include "G4Track.hh"
46 #include "G4ios.hh"
47 #include "G4SteppingManager.hh"
48 #include "G4ThreeVector.hh"
49
50 XrayTelAnalysis* XrayTelAnalysis::instance = 0
51
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() :
62 asciiFile(0),nEnteringTracks(0), totEntering
63 {
64 histFileName = "xraytel";
65
66
67 asciiFileName="xraytel.out";
68 asciiFile = new std::ofstream(asciiFileName)
69
70 if(asciiFile->is_open())
71 (*asciiFile) << "Energy (keV) x (mm) y
72 }
73
74 XrayTelAnalysis::~XrayTelAnalysis()
75 {
76 if (asciiFile)
77 delete asciiFile;
78 if (nEnteringTracks)
79 delete nEnteringTracks;
80 if (totEnteringEnergy)
81 delete totEnteringEnergy;
82 }
83
84
85 XrayTelAnalysis* XrayTelAnalysis::getInstance(
86 {
87 G4AutoLock l(&instanceMutex);
88 if (instance == 0)
89 instance = new XrayTelAnalysis;
90 return instance;
91 }
92
93
94 void XrayTelAnalysis::book(G4bool isMaster)
95 {
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
117 G4AnalysisManager* man = G4AnalysisManager::
118 man->SetDefaultFileType("root");
119
120 // Open an output file: it is done in master
121 // printout is done only by the master, for
122 if (isMaster)
123 G4cout << "Opening output file " << histFi
124 man->OpenFile(histFileName);
125 man->SetFirstHistoId(1);
126 if (isMaster)
127 G4cout << " done" << G4endl;
128
129 // Book 1D histograms
130 man->CreateH1("h1","Energy, all /keV", 100,
131 man->CreateH1("h2","Energy, entering detecto
132
133 // Book 2D histograms (notice: the numbering
134 man->CreateH2("d1","y-z, all /mm", 100,-500.
135 man->CreateH2("d2","y-z, entering detector /
136
137 // Book ntuples
138 man->CreateNtuple("tree", "Track ntuple");
139 man->CreateNtupleDColumn("energy");
140 man->CreateNtupleDColumn("x");
141 man->CreateNtupleDColumn("y");
142 man->CreateNtupleDColumn("z");
143 man->CreateNtupleDColumn("dirx");
144 man->CreateNtupleDColumn("diry");
145 man->CreateNtupleDColumn("dirz");
146 man->FinishNtuple();
147 }
148
149
150 void XrayTelAnalysis::finish(G4bool isMaster)
151 {
152 G4AutoLock l(&dataManipulationMutex);
153 // Save histograms
154 G4AnalysisManager* man = G4AnalysisManager::
155 man->Write();
156 man->CloseFile();
157 man->Clear();
158
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 }
231
232 void XrayTelAnalysis::analyseStepping(const G4
233 {
234 G4AutoLock l(&dataManipulationMutex);
235 eKin = track.GetKineticEnergy()/keV;
236 G4ThreeVector pos = track.GetPosition()/mm;
237 y = pos.y();
238 z = pos.z();
239 G4ThreeVector dir = track.GetMomentumDirecti
240 dirX = dir.x();
241 dirY = dir.y();
242 dirZ = dir.z();
243
244 // Fill histograms
245 G4AnalysisManager* man = G4AnalysisManager::
246 man->FillH1(1,eKin);
247 man->FillH2(1,y,z);
248
249 // Fill histograms and ntuple, tracks enteri
250 if (entering) {
251 // Fill and plot histograms
252 man->FillH1(2,eKin);
253 man->FillH2(2,y,z);
254
255 man->FillNtupleDColumn(0,eKin);
256 man->FillNtupleDColumn(1,x);
257 man->FillNtupleDColumn(2,y);
258 man->FillNtupleDColumn(3,z);
259 man->FillNtupleDColumn(4,dirX);
260 man->FillNtupleDColumn(5,dirY);
261 man->FillNtupleDColumn(6,dirZ);
262 man->AddNtupleRow();
263 }
264
265
266 // Write to file
267 if (entering) {
268 if(asciiFile->is_open()) {
269 (*asciiFile) << std::setiosflags(std::io
270 << std::setprecision(3)
271 << std::setiosflags(std::ios::right)
272 << std::setw(10);
273 (*asciiFile) << eKin;
274 (*asciiFile) << std::setiosflags(std::io
275 << std::setprecision(3)
276 << std::setiosflags(std::ios::right)
277 << std::setw(10);
278 (*asciiFile) << x;
279 (*asciiFile) << std::setiosflags(std::io
280 << std::setprecision(3)
281 << std::setiosflags(std::ios::right)
282 << std::setw(10);
283 (*asciiFile) << y;
284 (*asciiFile) << std::setiosflags(std::io
285 << std::setprecision(3)
286 << std::setiosflags(std::ios::right)
287 << std::setw(10);
288 (*asciiFile) << z
289 << G4endl;
290 }
291 }
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 }
316
317