source: trunk/examples/advanced/xray_telescope/src/XrayTelAnalysis.cc@ 1327

Last change on this file since 1327 was 1230, checked in by garnier, 16 years ago

update to geant4.9.3
File size: 8.8 KB
RevLine 
[807]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// $Id: XrayTelAnalysis.cc,v 1.12 2006/06/29 16:30:00 gunter Exp $
[1230]28// GEANT4 tag $Name: geant4-09-03-cand-01 $
[807]29//
30// Author: A. Pfeiffer (Andreas.Pfeiffer@cern.ch)
31// (copied from his UserAnalyser class)
32//
33// History:
34// -----------
35// 8 Nov 2002 GS Migration to AIDA 3
36// 7 Nov 2001 MGP Implementation
37//
38// -------------------------------------------------------------------
39
40#include "XrayTelAnalysis.hh"
41#include "globals.hh"
42#include "G4Track.hh"
43#include "G4ios.hh"
44#include <fstream>
45#include <iomanip>
46#include "G4SteppingManager.hh"
47#include "G4ThreeVector.hh"
48
49
50XrayTelAnalysis* XrayTelAnalysis::instance = 0;
51
52XrayTelAnalysis::XrayTelAnalysis()
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{
64#ifdef G4ANALYSIS_USE
65 histFileName = "xraytel";
66 histFileType = "hbook";
67#endif
68
69 asciiFileName="xraytel.out";
70 std::ofstream asciiFile(asciiFileName, std::ios::app);
71 if(asciiFile.is_open()) {
72 asciiFile << "Energy (keV) x (mm) y (mm) z (mm)" << G4endl << G4endl;
73 }
74}
75
76XrayTelAnalysis::~XrayTelAnalysis()
77{
78#ifdef G4ANALYSIS_USE
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}
98
99XrayTelAnalysis* XrayTelAnalysis::getInstance()
100{
101 if (instance == 0) instance = new XrayTelAnalysis;
102 return instance;
103}
104
105
106void XrayTelAnalysis::book()
107{
108#ifdef G4ANALYSIS_USE
109 //build up the factories
110 analysisFactory = AIDA_createAnalysisFactory();
111 if(analysisFactory) {
112 //parameters for the TreeFactory
113 G4bool fileExists = true;
114 G4bool readOnly = false;
115 AIDA::ITreeFactory* treeFactory = analysisFactory->createTreeFactory();
116 if(treeFactory) {
117 G4String histFileNameComplete;
118 histFileNameComplete = histFileName+".hbook";
119 tree = treeFactory->create(histFileNameComplete, "hbook", readOnly, fileExists);
120 G4cout << " Histogramfile: " << histFileNameComplete << G4endl;
121
122 if (tree) {
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" );
139 }
140 delete treeFactory;
141 }
142 }
143#endif
144
145}
146
147void XrayTelAnalysis::finish()
148{
149#ifdef G4ANALYSIS_USE
150 if (tree) {
151 // Committing the transaction with the tree
152 G4cout << "Committing..." << G4endl;
153 // write all histograms to file
154 tree->commit();
155
156 G4cout << "Closing the tree..." << G4endl;
157
158 // close (will again commit)
159 tree->close();
160 }
161
162 // extra delete as objects are created in book() method rather than during
163 // initialisation of class
164
165// #ifdef G4ANALYSIS_USE_PLOTTER
166// if (plotterFactory) delete plotterFactory;
167// #endif
168
169 if (tupleFactory) delete tupleFactory;
170 if (histoFactory) delete histoFactory;
171 if (tree) delete tree;
172 if (analysisFactory) delete analysisFactory;
173#endif
174}
175
176void XrayTelAnalysis::analyseStepping(const G4Track& track, G4bool entering)
177{
178 eKin = track.GetKineticEnergy()/keV;
179 G4ThreeVector pos = track.GetPosition()/mm;
180 y = pos.y();
181 z = pos.z();
182 G4ThreeVector dir = track.GetMomentumDirection();
183 dirX = dir.x();
184 dirY = dir.y();
185 dirZ = dir.z();
186
187#ifdef G4ANALYSIS_USE
188 // Fill histograms, all tracks
189
190 h1->fill(eKin); // fill(x,y,weight)
191
192 h2->fill(y,z);
193
194 // Fill histograms and ntuple, tracks entering the detector
195 if (entering) {
196 // Fill and plot histograms
197
198 h3->fill(eKin);
199
200 h4->fill(y,z);
201// #ifdef G4ANALYSIS_USE_PLOTTER
202// plotAll();
203// #endif
204 }
205
206 // Fill ntuple
207 if (entering) {
208 if (ntuple) {
209 // Fill the secondaries ntuple
210 ntuple->fill( ntuple->findColumn( "energy" ), (G4double) eKin );
211 ntuple->fill( ntuple->findColumn( "x" ), (G4double) x );
212 ntuple->fill( ntuple->findColumn( "y" ), (G4double) y );
213 ntuple->fill( ntuple->findColumn( "z" ), (G4double) z );
214 ntuple->fill( ntuple->findColumn( "dirx" ), (G4double) dirX );
215 ntuple->fill( ntuple->findColumn( "diry" ), (G4double) dirY );
216 ntuple->fill( ntuple->findColumn( "dirz" ), (G4double) dirZ );
217
218 ntuple->addRow(); // check for returning true ...
219 } else {
220 G4cout << "Ntuple not found" << G4endl;
221 }
222 }
223
224#endif
225
226 // Write to file
227 if (entering) {
228 std::ofstream asciiFile(asciiFileName, std::ios::app);
229 if(asciiFile.is_open()) {
230 asciiFile << std::setiosflags(std::ios::fixed)
231 << std::setprecision(3)
232 << std::setiosflags(std::ios::right)
233 << std::setw(10);
234 asciiFile << eKin;
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();
252 }
253 }
254
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
293
Note: See TracBrowser for help on using the repository browser.