| 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,v 1.12 2006/06/29 16:30:00 gunter Exp $
>> 28 // GEANT4 tag $Name: geant4-09-01 $
27 // 29 //
28 // Author: A. Pfeiffer (Andreas.Pfeiffer@cern 30 // Author: A. Pfeiffer (Andreas.Pfeiffer@cern.ch)
29 // (copied from his UserAnalyser class 31 // (copied from his UserAnalyser class)
30 // 32 //
31 // History: 33 // History:
32 // ----------- 34 // -----------
33 // 19 Mar 2013 LP Migrated to G4tool <<
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 <iomanip> <<
41 #include "G4AutoLock.hh" <<
42 #include "XrayTelAnalysis.hh" 40 #include "XrayTelAnalysis.hh"
43 #include "globals.hh" 41 #include "globals.hh"
44 #include "G4SystemOfUnits.hh" <<
45 #include "G4Track.hh" 42 #include "G4Track.hh"
46 #include "G4ios.hh" 43 #include "G4ios.hh"
>> 44 #include <fstream>
>> 45 #include <iomanip>
47 #include "G4SteppingManager.hh" 46 #include "G4SteppingManager.hh"
48 #include "G4ThreeVector.hh" 47 #include "G4ThreeVector.hh"
49 48
50 XrayTelAnalysis* XrayTelAnalysis::instance = 0 <<
51 49
52 namespace { << 50 XrayTelAnalysis* XrayTelAnalysis::instance = 0;
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 51
61 XrayTelAnalysis::XrayTelAnalysis() : << 52 XrayTelAnalysis::XrayTelAnalysis()
62 asciiFile(0),nEnteringTracks(0), totEntering << 53 #ifdef G4ANALYSIS_USE
>> 54 : analysisFactory(0)
>> 55 , tree(0)
>> 56 , histoFactory(0)
>> 57 , tupleFactory(0), h1(0), h2(0), h3(0), h4(0), ntuple(0)
>> 58 // #ifdef G4ANALYSIS_USE_PLOTTER
>> 59 // , plotterFactory(0)
>> 60 // , plotter(0)
>> 61 // #endif
>> 62 #endif
63 { 63 {
>> 64 #ifdef G4ANALYSIS_USE
64 histFileName = "xraytel"; 65 histFileName = "xraytel";
65 << 66 histFileType = "hbook";
>> 67 #endif
66 68
67 asciiFileName="xraytel.out"; 69 asciiFileName="xraytel.out";
68 asciiFile = new std::ofstream(asciiFileName) << 70 std::ofstream asciiFile(asciiFileName, std::ios::app);
69 << 71 if(asciiFile.is_open()) {
70 if(asciiFile->is_open()) << 72 asciiFile << "Energy (keV) x (mm) y (mm) z (mm)" << G4endl << G4endl;
71 (*asciiFile) << "Energy (keV) x (mm) y << 73 }
72 } 74 }
73 75
74 XrayTelAnalysis::~XrayTelAnalysis() 76 XrayTelAnalysis::~XrayTelAnalysis()
75 { << 77 {
76 if (asciiFile) << 78 #ifdef G4ANALYSIS_USE
77 delete asciiFile; <<
78 if (nEnteringTracks) <<
79 delete nEnteringTracks; <<
80 if (totEnteringEnergy) <<
81 delete totEnteringEnergy; <<
82 } <<
83 79
>> 80 // #ifdef G4ANALYSIS_USE_PLOTTER
>> 81 // if (plotterFactory) delete plotterFactory;
>> 82 // plotterFactory = 0;
>> 83 // #endif
>> 84
>> 85 if (tupleFactory) delete tupleFactory;
>> 86 tupleFactory = 0;
>> 87
>> 88 if (histoFactory) delete histoFactory;
>> 89 histoFactory = 0;
>> 90
>> 91 if (tree) delete tree;
>> 92 tree = 0;
>> 93
>> 94 if (analysisFactory) delete analysisFactory;
>> 95 analysisFactory = 0;
>> 96 #endif
>> 97 }
84 98
85 XrayTelAnalysis* XrayTelAnalysis::getInstance( 99 XrayTelAnalysis* XrayTelAnalysis::getInstance()
86 { 100 {
87 G4AutoLock l(&instanceMutex); << 101 if (instance == 0) instance = new XrayTelAnalysis;
88 if (instance == 0) <<
89 instance = new XrayTelAnalysis; <<
90 return instance; 102 return instance;
91 } 103 }
92 104
93 105
94 void XrayTelAnalysis::book(G4bool isMaster) << 106 void XrayTelAnalysis::book()
95 { 107 {
96 G4AutoLock l(&dataManipulationMutex); << 108 #ifdef G4ANALYSIS_USE
97 << 109 //build up the factories
98 //reset counters: do be done only once, by t << 110 analysisFactory = AIDA_createAnalysisFactory();
99 if (isMaster) << 111 if(analysisFactory) {
100 { << 112 //parameters for the TreeFactory
101 if (nEnteringTracks) << 113 G4bool fileExists = true;
102 { << 114 G4bool readOnly = false;
103 delete nEnteringTracks; << 115 AIDA::ITreeFactory* treeFactory = analysisFactory->createTreeFactory();
104 nEnteringTracks = 0; << 116 if(treeFactory) {
105 } << 117 G4String histFileNameComplete;
106 nEnteringTracks = new std::map<G4int,G4i << 118 histFileNameComplete = histFileName+".hbook";
107 << 119 tree = treeFactory->create(histFileNameComplete, "hbook", readOnly, fileExists);
108 if (totEnteringEnergy) << 120 G4cout << " Histogramfile: " << histFileNameComplete << G4endl;
109 { << 121
110 delete totEnteringEnergy; << 122 if (tree) {
111 totEnteringEnergy = 0; << 123 G4cout << "Tree store : " << tree->storeName() << G4endl;
>> 124 G4cout << "Booked Hbook File " << G4endl;
>> 125
>> 126 //HistoFactory and TupleFactory depend on theTree
>> 127 histoFactory = analysisFactory->createHistogramFactory ( *tree );
>> 128 tupleFactory = analysisFactory->createTupleFactory ( *tree );
>> 129
>> 130 // Book histograms
>> 131 h1 = histoFactory->createHistogram1D("1","Energy, all /keV", 100,0.,100.);
>> 132 h2 = histoFactory->createHistogram2D("2","y-z, all /mm", 100,-500.,500.,100,-500.,500.);
>> 133 h3 = histoFactory->createHistogram1D("3","Energy, entering detector /keV", 500,0.,500.);
>> 134 h4 = histoFactory->createHistogram2D("4","y-z, entering detector /mm", 200,-50.,50.,200,-50.,50.);
>> 135
>> 136 // Book ntuples
>> 137 ntuple = tupleFactory->create( "10", "Track ntuple",
>> 138 "double energy,x,y,z,dirx,diry,dirz" );
112 } 139 }
113 totEnteringEnergy = new std::map<G4int,G << 140 delete treeFactory;
114 } 141 }
>> 142 }
>> 143 #endif
115 144
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 } 145 }
148 146
149 << 147 void XrayTelAnalysis::finish()
150 void XrayTelAnalysis::finish(G4bool isMaster) <<
151 { 148 {
152 G4AutoLock l(&dataManipulationMutex); << 149 #ifdef G4ANALYSIS_USE
153 // Save histograms << 150 if (tree) {
154 G4AnalysisManager* man = G4AnalysisManager:: << 151 // Committing the transaction with the tree
155 man->Write(); << 152 G4cout << "Committing..." << G4endl;
156 man->CloseFile(); << 153 // write all histograms to file
157 man->Clear(); << 154 tree->commit();
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 155
179 //MT run: sum results << 156 G4cout << "Closing the tree..." << G4endl;
180 <<
181 157
182 //MT build, but sequential run << 158 // close (will again commit)
183 if (nEnteringTracks->count(-1)) << 159 tree->close();
184 { << 160 }
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 161
196 G4bool loopAgain = true; << 162 // extra delete as objects are created in book() method rather than during
197 G4int totEntries = 0; << 163 // initialisation of class
198 G4double totEnergy = 0.; << 164
199 << 165 // #ifdef G4ANALYSIS_USE_PLOTTER
200 G4cout << "################################# << 166 // if (plotterFactory) delete plotterFactory;
201 for (size_t i=0; loopAgain; i++) << 167 // #endif
202 { << 168
203 //ok, this thread was found << 169 if (tupleFactory) delete tupleFactory;
204 if (nEnteringTracks->count(i)) << 170 if (histoFactory) delete histoFactory;
205 { << 171 if (tree) delete tree;
206 G4cout << "End of Run summary (thread= " < << 172 if (analysisFactory) delete analysisFactory;
207 G4int part = nEnteringTracks->find(i)->sec << 173 #endif
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 } 174 }
231 175
232 void XrayTelAnalysis::analyseStepping(const G4 176 void XrayTelAnalysis::analyseStepping(const G4Track& track, G4bool entering)
233 { 177 {
234 G4AutoLock l(&dataManipulationMutex); <<
235 eKin = track.GetKineticEnergy()/keV; 178 eKin = track.GetKineticEnergy()/keV;
236 G4ThreeVector pos = track.GetPosition()/mm; 179 G4ThreeVector pos = track.GetPosition()/mm;
237 y = pos.y(); 180 y = pos.y();
238 z = pos.z(); 181 z = pos.z();
239 G4ThreeVector dir = track.GetMomentumDirecti 182 G4ThreeVector dir = track.GetMomentumDirection();
240 dirX = dir.x(); 183 dirX = dir.x();
241 dirY = dir.y(); 184 dirY = dir.y();
242 dirZ = dir.z(); 185 dirZ = dir.z();
243 186
244 // Fill histograms << 187 #ifdef G4ANALYSIS_USE
245 G4AnalysisManager* man = G4AnalysisManager:: << 188 // Fill histograms, all tracks
246 man->FillH1(1,eKin); << 189
247 man->FillH2(1,y,z); << 190 h1->fill(eKin); // fill(x,y,weight)
248 << 191
>> 192 h2->fill(y,z);
>> 193
249 // Fill histograms and ntuple, tracks enteri 194 // Fill histograms and ntuple, tracks entering the detector
250 if (entering) { 195 if (entering) {
251 // Fill and plot histograms 196 // Fill and plot histograms
252 man->FillH1(2,eKin); <<
253 man->FillH2(2,y,z); <<
254 197
255 man->FillNtupleDColumn(0,eKin); << 198 h3->fill(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 199
>> 200 h4->fill(y,z);
>> 201 // #ifdef G4ANALYSIS_USE_PLOTTER
>> 202 // plotAll();
>> 203 // #endif
>> 204 }
265 205
266 // Write to file << 206 // Fill ntuple
267 if (entering) { 207 if (entering) {
268 if(asciiFile->is_open()) { << 208 if (ntuple) {
269 (*asciiFile) << std::setiosflags(std::io << 209 // Fill the secondaries ntuple
270 << std::setprecision(3) << 210 ntuple->fill( ntuple->findColumn( "energy" ), (G4double) eKin );
271 << std::setiosflags(std::ios::right) << 211 ntuple->fill( ntuple->findColumn( "x" ), (G4double) x );
272 << std::setw(10); << 212 ntuple->fill( ntuple->findColumn( "y" ), (G4double) y );
273 (*asciiFile) << eKin; << 213 ntuple->fill( ntuple->findColumn( "z" ), (G4double) z );
274 (*asciiFile) << std::setiosflags(std::io << 214 ntuple->fill( ntuple->findColumn( "dirx" ), (G4double) dirX );
275 << std::setprecision(3) << 215 ntuple->fill( ntuple->findColumn( "diry" ), (G4double) dirY );
276 << std::setiosflags(std::ios::right) << 216 ntuple->fill( ntuple->findColumn( "dirz" ), (G4double) dirZ );
277 << std::setw(10); << 217
278 (*asciiFile) << x; << 218 ntuple->addRow(); // check for returning true ...
279 (*asciiFile) << std::setiosflags(std::io << 219 } else {
280 << std::setprecision(3) << 220 G4cout << "Ntuple not found" << G4endl;
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 } 221 }
291 } 222 }
292 } <<
293 223
294 void XrayTelAnalysis::Update(G4double energy,G << 224 #endif
295 { << 225
296 G4AutoLock l(&dataManipulationMutex); << 226 // Write to file
297 //It already exists: increase the counter << 227 if (entering) {
298 if (nEnteringTracks->count(threadID)) << 228 std::ofstream asciiFile(asciiFileName, std::ios::app);
299 { << 229 if(asciiFile.is_open()) {
300 (nEnteringTracks->find(threadID)->second << 230 asciiFile << std::setiosflags(std::ios::fixed)
301 } << 231 << std::setprecision(3)
302 else //enter a new one << 232 << std::setiosflags(std::ios::right)
303 { << 233 << std::setw(10);
304 G4int tracks = 1; << 234 asciiFile << eKin;
305 nEnteringTracks->insert(std::make_pair(t << 235 asciiFile << std::setiosflags(std::ios::fixed)
>> 236 << std::setprecision(3)
>> 237 << std::setiosflags(std::ios::right)
>> 238 << std::setw(10);
>> 239 asciiFile << x;
>> 240 asciiFile << std::setiosflags(std::ios::fixed)
>> 241 << std::setprecision(3)
>> 242 << std::setiosflags(std::ios::right)
>> 243 << std::setw(10);
>> 244 asciiFile << y;
>> 245 asciiFile << std::setiosflags(std::ios::fixed)
>> 246 << std::setprecision(3)
>> 247 << std::setiosflags(std::ios::right)
>> 248 << std::setw(10);
>> 249 asciiFile << z
>> 250 << G4endl;
>> 251 asciiFile.close();
306 } 252 }
>> 253 }
307 254
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 } 255 }
>> 256
>> 257 // #ifdef G4ANALYSIS_USE_PLOTTER
>> 258 // void XrayTelAnalysis::plotAll()
>> 259 // {
>> 260 // if (!plotter) {
>> 261 // AIDA::IPlotterFactory* plotterFactory =
>> 262 // analysisFactory->createPlotterFactory();
>> 263 // if(plotterFactory) {
>> 264 // G4cout << "Creating the Plotter" << G4endl;
>> 265 // plotter = plotterFactory->create();
>> 266 // if(plotter) {
>> 267 // // Map the plotter on screen :
>> 268 // G4cout << "Showing the Plotter on screen" << G4endl;
>> 269 // plotter->show();
>> 270 // } else {
>> 271 // G4cout << "XrayTelAnalysis::plotAll: WARNING: Plotter not created" << G4endl;
>> 272 // }
>> 273 // delete plotterFactory;
>> 274 // } else {
>> 275 // G4cout << "XrayTelAnalysis::plotAll: WARNING: Plotter Factory not created" << G4endl;
>> 276 // }
>> 277 // }
>> 278
>> 279 // if (plotter) {
>> 280 // plotter->createRegions(2,1,0);
>> 281 // AIDA::IHistogram1D* hp = dynamic_cast<AIDA::IHistogram1D *> ( tree->find("3") );
>> 282 // AIDA::IHistogram1D& h = *hp;
>> 283 // (plotter->currentRegion()).plot(h);
>> 284 // plotter->refresh();
>> 285 // plotter->setCurrentRegionNumber(1);
>> 286 // AIDA::IHistogram1D* hp2 = dynamic_cast<AIDA::IHistogram1D *> ( tree->find("1") );
>> 287 // AIDA::IHistogram1D& h2 = *hp2;
>> 288 // (plotter->currentRegion()).plot(h2);
>> 289 // plotter->refresh();
>> 290 // }
>> 291 // }
>> 292 // #endif
316 293
317 294