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