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

Last change on this file since 1358 was 1337, checked in by garnier, 14 years ago

tag geant4.9.4 beta 1 + modifs locales
File size: 8.8 KB
Line 
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 $
28// GEANT4 tag $Name: geant4-09-04-beta-01 $
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.