source: trunk/examples/advanced/gammaray_telescope/src/GammaRayTelAnalysis.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: 11.5 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#ifdef G4ANALYSIS_USE
27//
28// $Id: GammaRayTelAnalysis.cc,v 1.21 2006/06/29 15:56:07 gunter Exp $
29// GEANT4 tag $Name: geant4-09-03-cand-01 $
30// ------------------------------------------------------------
31// GEANT 4 class implementation file
32// CERN Geneva Switzerland
33//
34//
35// ------------ GammaRayAnalysisManager ------
36// by R.Giannitrapani, F.Longo & G.Santin (03 dic 2000)
37//
38// 29.05.2003 F.Longo
39// - anaphe 5.0.5 compliant
40//
41// 18.06.2002 R.Giannitrapani, F.Longo & G.Santin
42// - new release for Anaphe 4.0.3
43//
44// 07.12.2001 A.Pfeiffer
45// - integrated Guy's addition of the ntuple
46//
47// 06.12.2001 A.Pfeiffer
48// - updating to new design (singleton)
49//
50// 22.11.2001 G.Barrand
51// - Adaptation to AIDA
52//
53// ************************************************************
54#include <fstream>
55
56#include "G4RunManager.hh"
57
58#include "GammaRayTelAnalysis.hh"
59#include "GammaRayTelDetectorConstruction.hh"
60#include "GammaRayTelAnalysisMessenger.hh"
61
62GammaRayTelAnalysis* GammaRayTelAnalysis::instance = 0;
63
64//--------------------------------------------------------------------------------
65//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66
67GammaRayTelAnalysis::GammaRayTelAnalysis()
68 :GammaRayTelDetector(0),analysisFactory(0), tree(0)//, plotter(0),
69 ,tuple(0)
70 ,energy(0), hits(0), posXZ(0), posYZ(0)
71 ,histo1DDraw("enable"),histo1DSave("enable"),histo2DDraw("enable")
72 ,histo2DSave("enable"),histo2DMode("strip")
73{
74 G4RunManager* runManager = G4RunManager::GetRunManager();
75 GammaRayTelDetector =
76 (GammaRayTelDetectorConstruction*)(runManager->GetUserDetectorConstruction());
77
78#ifdef G4ANALYSIS_USE
79 // Define the messenger and the analysis system
80
81 analysisMessenger = new GammaRayTelAnalysisMessenger(this);
82
83 analysisFactory = AIDA_createAnalysisFactory(); // create the Analysis Factory
84 if(analysisFactory) {
85 AIDA::ITreeFactory* treeFactory = analysisFactory->createTreeFactory();
86 // create Tree Factory
87
88 if(treeFactory) {
89 // Tree in memory :
90 // Create a "tree" associated to an xml file
91
92 // tree = treeFactory->create("gammaraytel.hbook", "hbook", false, false);
93 // (hbook implementation)
94 tree = treeFactory->create("gammaraytel.aida","xml",false,true,"uncompress");
95 if(tree) {
96 // Get a tuple factory :
97 ITupleFactory* tupleFactory = analysisFactory->createTupleFactory(*tree);
98 if(tupleFactory) {
99 // Create a tuple :
100 tuple = tupleFactory->create("1","1", "float energy, plane, x, y, z");
101 assert(tuple);
102
103 delete tupleFactory;
104 }
105
106 IHistogramFactory* histoFactory = analysisFactory->createHistogramFactory(*tree);
107
108 if(histoFactory) {
109 // Create histos :
110
111 int Nplane = GammaRayTelDetector->GetNbOfTKRLayers();
112 int Nstrip = GammaRayTelDetector->GetNbOfTKRStrips();
113 int Ntile = GammaRayTelDetector->GetNbOfTKRTiles();
114 double sizexy = GammaRayTelDetector->GetTKRSizeXY();
115 double sizez = GammaRayTelDetector->GetTKRSizeZ();
116 int N = Nstrip*Ntile;
117
118 // 1D histogram that store the energy deposition of the
119 // particle in the last (number 0) TKR X-plane
120
121 energy = histoFactory->createHistogram1D("10","Edep in the last X plane (keV)", 100, 50, 200);
122
123 // 1D histogram that store the hits distribution along the TKR X-planes
124
125 hits = histoFactory->createHistogram1D("20","Hits dist in TKR X planes",Nplane, 0, Nplane-1);
126 // 2D histogram that store the position (mm) of the hits (XZ projection)
127
128 if (histo2DMode == "strip"){
129 posXZ = histoFactory->createHistogram2D("30","Tracker Hits XZ (strip,plane)",
130 N, 0, N-1,
131 2*Nplane, 0, Nplane-1);
132 }
133 else
134 {
135 posXZ = histoFactory->createHistogram2D("30","Tracker Hits XZ (x,z) in mm",
136 int(sizexy/5), -sizexy/2, sizexy/2,
137 int(sizez/5), -sizez/2, sizez/2);
138 }
139 // 2D histogram that store the position (mm) of the hits (YZ projection)
140
141 if(histo2DMode=="strip")
142 posYZ = histoFactory->createHistogram2D("40","Tracker Hits YZ (strip,plane)",
143 N, 0, N-1,
144 2*Nplane, 0, Nplane-1);
145 else
146 posYZ = histoFactory->createHistogram2D("40","Tracker Hits YZ (y,z) in mm",
147 int(sizexy/5), -sizexy/2, sizexy/2,
148 int(sizez/5), -sizez/2, sizez/2);
149
150 delete histoFactory;
151 }
152
153 }
154 }
155 delete treeFactory; // Will not delete the ITree.
156
157 }
158
159 // IPlotterFactory* plotterFactory = analysisFactory->createPlotterFactory(0,0);
160// if(plotterFactory) {
161// plotter = plotterFactory->create();
162// if(plotter) {
163// plotter->show();
164// plotter->setParameter("pageTitle","Gamma Ray Tel");
165
166// }
167// delete plotterFactory;
168// }
169
170#endif
171}
172
173
174GammaRayTelAnalysis::~GammaRayTelAnalysis() {
175 Finish();
176}
177
178
179void GammaRayTelAnalysis::Init()
180{
181}
182
183void GammaRayTelAnalysis::Finish()
184{
185#ifdef G4ANALYSIS_USE
186 delete tree;
187 //delete plotter;
188 // delete analysisFactory; // Will delete tree and histos.
189 delete analysisMessenger;
190
191 analysisMessenger = 0;
192#endif
193}
194
195GammaRayTelAnalysis* GammaRayTelAnalysis::getInstance()
196{
197 if (instance == 0) instance = new GammaRayTelAnalysis();
198 return instance;
199}
200
201// This function fill the 2d histogram of the XZ positions
202void GammaRayTelAnalysis::InsertPositionXZ(double x, double z)
203{
204#ifdef G4ANALYSIS_USE
205 if(posXZ) posXZ->fill(x, z);
206#endif
207}
208
209// This function fill the 2d histogram of the YZ positions
210void GammaRayTelAnalysis::InsertPositionYZ(double y, double z)
211{
212#ifdef G4ANALYSIS_USE
213 if(posYZ) posYZ->fill(y, z);
214#endif
215}
216
217// This function fill the 1d histogram of the energy released in the last Si plane
218void GammaRayTelAnalysis::InsertEnergy(double en)
219{
220#ifdef G4ANALYSIS_USE
221 if(energy) energy->fill(en);
222#endif
223}
224
225// This function fill the 1d histogram of the hits distribution along the TKR planes
226void GammaRayTelAnalysis::InsertHits(int nplane)
227{
228#ifdef G4ANALYSIS_USE
229 if(hits) hits->fill(nplane);
230#endif
231}
232
233void GammaRayTelAnalysis::setNtuple(float E, float p, float x, float y, float z)
234{
235#ifdef G4ANALYSIS_USE
236 tuple->fill(tuple->findColumn("energy"),E);
237 tuple->fill(tuple->findColumn("plane"),p);
238 tuple->fill(tuple->findColumn("x"),x);
239 tuple->fill(tuple->findColumn("y"),y);
240 tuple->fill(tuple->findColumn("z"),z);
241 tuple->addRow();
242#endif
243}
244
245//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
246/*
247 This member reset the histograms and it is called at the begin
248 of each run; here we put the inizialization so that the histograms have
249 always the right dimensions depending from the detector geometry
250*/
251
252//void GammaRayTelAnalysis::BeginOfRun(G4int n) // to be reintroduced
253void GammaRayTelAnalysis::BeginOfRun()
254{
255#ifdef G4ANALYSIS_USE
256
257 if(energy) energy->reset();
258 if(hits) hits->reset();
259 if(posXZ) posXZ->reset();
260 if(posYZ) posYZ->reset();
261
262#endif
263}
264
265//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
266
267/*
268 This member is called at the end of each run
269*/
270void GammaRayTelAnalysis::EndOfRun(G4int)
271{
272#ifdef G4ANALYSIS_USE
273 if(tree) {
274 tree->commit();
275 tree->close();
276 }
277
278// if(plotter) {
279
280// // We set one single region for the plotter
281// // We now print the histograms, each one in a separate file
282
283// if(histo2DSave == "enable") {
284// char name[15];
285// plotter->createRegions(1,1);
286// sprintf(name,"posxz_%d.ps", n);
287// plotter->currentRegion().plot(*posXZ);
288// plotter->refresh();
289// // plotter->write(name,"ps"); // temporary unavailable
290
291// plotter->createRegions(1,1);
292// sprintf(name,"posyz_%d.ps", n);
293// plotter->currentRegion().plot(*posYZ);
294// plotter->next().plot(*posYZ);
295// plotter->refresh();
296// // plotter->write(name,"ps"); // temporary unavailable
297// }
298
299// if(histo1DSave == "enable") {
300// plotter->createRegions(1,1);
301// char name[15];
302// sprintf(name,"energy_%d.ps", n);
303// plotter->currentRegion().plot(*energy);
304// plotter->refresh();
305// // plotter->write(name,"ps"); // temporary unavailable
306// plotter->createRegions(1,1);
307// sprintf(name,"hits_%d.ps", n);
308// plotter->currentRegion().plot(*hits);
309// plotter->refresh();
310// // plotter->write(name,"ps"); // temporary unavailable
311// plotter->createRegions(1,2);
312// plotter->currentRegion().plot(*energy);
313// plotter->next().plot(*hits);
314// plotter->refresh();
315 // }
316#endif
317}
318
319/* This member is called at the end of every event */
320void GammaRayTelAnalysis::EndOfEvent(G4int flag)
321{
322 // The plotter is updated only if there is some
323 // hits in the event
324 if(!flag) return;
325#ifdef G4ANALYSIS_USE
326 // Set the plotter ; set the number of regions and attach histograms
327 // to plot for each region.
328 // It is done here, since then EndOfRun set regions
329 // for paper output.
330 // if(plotter) {
331// if((histo2DDraw == "enable") && (histo1DDraw == "enable")) {
332// plotter->createRegions(1,2);
333// //plotter->currentRegion().plot(*posXZ); //temporary unavailable
334// plotter->currentRegion().plot(*hits);
335// // plotter->next().plot(*posYZ); //temporary unavailable
336// plotter->next().plot(*energy);
337// //plotter->next().plot(*energy);
338// // plotter->currentRegion().plot(*hits);
339// //plotter->next().plot(*hits);
340// } else if((histo1DDraw == "enable") && (histo2DDraw != "enable")) {
341// plotter->createRegions(1,2);
342// plotter->currentRegion().plot(*energy);
343// plotter->next().plot(*hits);
344// } else if((histo1DDraw != "enable") && (histo2DDraw == "enable")) {
345// /* plotter->createRegions(1,2);
346// plotter->currentRegion().plot(*posXZ);
347// plotter->next().plot(*posYZ);*/
348// G4cout << "Temporary Unavailable " << G4endl;
349// } else { // Nothing to plot.
350// plotter->createRegions(1,1);
351// }
352// plotter->refresh();
353// }
354
355#endif
356}
357#endif
358
359
360
361
362
363
364
Note: See TracBrowser for help on using the repository browser.