Geant4_10
HistoManager.cc
Go to the documentation of this file.
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 // $Id: HistoManager.cc 73008 2013-08-15 08:43:13Z gcosmo $
27 //
28 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
29 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
30 
31 #include "HistoManager.hh"
32 #include "HistoMessenger.hh"
33 #include "G4UnitsTable.hh"
34 
35 #ifdef G4ANALYSIS_USE
36 #include "AIDA/AIDA.h"
37 #endif
38 
39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
40 
42 #ifdef G4ANALYSIS_USE
43 : factoryOn(false),af(0),tree(0)
44 #endif
45 {
46 #ifdef G4ANALYSIS_USE
47  // Creating the analysis factory
48  af = AIDA_createAnalysisFactory();
49  if(!af) {
50  G4cout << " HistoManager::HistoManager :"
51  << " problem creating the AIDA analysis factory."
52  << G4endl;
53  }
54 #endif
55 
56  fileName[0] = "ecal";
57  fileType = "root";
58  fileOption = "export=root";
59 
60  // histograms
61  for (G4int k=0; k<MaxHisto; k++) {
62  histo[k] = 0;
63  exist[k] = false;
64  Unit[k] = 1.0;
65  Width[k] = 1.0;
66  ascii[k] = false;
67  }
68 
69  // ntuples
70  for (G4int k=0; k<MaxNtupl; k++) {
71  ntupl[k] = 0;
72  existNt[k] = false;
73  }
74 
75  histoMessenger = new HistoMessenger(this);
76 }
77 
78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
79 
81 {
82  delete histoMessenger;
83 
84 #ifdef G4ANALYSIS_USE
85  delete af;
86 #endif
87 }
88 
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90 
92 {
93 #ifdef G4ANALYSIS_USE
94  if(!af) return;
95 
96  // Creating a tree mapped to an histogram file.
97  fileName[1] = fileName[0] + "." + fileType;
98  G4bool readOnly = false;
99  G4bool createNew = true;
100  AIDA::ITreeFactory* tf = af->createTreeFactory();
101  tree = tf->create(fileName[1], fileType, readOnly, createNew, fileOption);
102  delete tf;
103  if(!tree) {
104  G4cout << " HistoManager::book :"
105  << " problem creating the AIDA tree with "
106  << " storeName = " << fileName[1]
107  << " storeType = " << fileType
108  << " readOnly = " << readOnly
109  << " createNew = " << createNew
110  << " options = " << fileOption
111  << G4endl;
112  return;
113  }
114 
115  // Creating a histogram factory, whose histograms will be handled by the tree
116  AIDA::IHistogramFactory* hf = af->createHistogramFactory(*tree);
117 
118  // create selected histograms
119  for (G4int k=1; k<MaxHisto; k++) {
120  if (exist[k]) {
121  histo[k] = hf->createHistogram1D( Label[k], Title[k],
122  Nbins[k], Vmin[k], Vmax[k]);
123  factoryOn = true;
124  }
125  }
126  delete hf;
127 
128  // Creating a ntuple factory, handled by the tree
129  AIDA::ITupleFactory* ntf = af->createTupleFactory(*tree);
130 
131  // create selected ntuples
132  for (G4int k=1; k<MaxNtupl; k++) {
133  if (existNt[k]) {
134  ntupl[k] = ntf->create( LabelNt[k], TitleNt[k], ColumnNt[k]);
135  factoryOn = true;
136  }
137  }
138 
139  delete ntf;
140 
141  if (factoryOn)
142  G4cout << "\n----> Histogram Tree is opened in " << fileName[1] << G4endl;
143 #endif
144 }
145 
146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
147 
149 {
150 #ifdef G4ANALYSIS_USE
151  if (factoryOn) {
152  saveAscii(); // Write ascii file, if any
153  tree->commit(); // Writing the histograms to the file
154  tree->close(); // and closing the tree (and the file)
155  G4cout << "\n----> Histogram Tree is saved in " << fileName[1] << G4endl;
156 
157  delete tree;
158  tree = 0;
159  factoryOn = false;
160  }
161 #endif
162 }
163 
164 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
165 
167 {
168  if (ih >= MaxHisto) {
169  G4cout << "---> warning from HistoManager::FillHisto() : histo " << ih
170  << " xbin= " << xbin << " weight= " << weight << G4endl;
171  return;
172  }
173 #ifdef G4ANALYSIS_USE
174  if(exist[ih]) histo[ih]->fill(xbin/Unit[ih], weight);
175 #endif
176 }
177 
178 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
179 
181  G4int nbins, G4double valmin, G4double valmax, const G4String& unit)
182 {
183  if (ih < 1 || ih >= MaxHisto) {
184  G4cout << "---> warning from HistoManager::SetHisto() : histo " << ih
185  << "does not exist" << G4endl;
186  return;
187  }
188 
189  const G4String id[] = { "0", "1", "2", "3", "4", "5" };
190  const G4String title[] =
191  { "dummy", //0
192  "total Evis in Ecal", //1
193  "total Edep in Ecal", //2
194  "Evis profile", //3
195  "Edep profile", //4
196  "Nb of Radiation Length" //5
197  };
198 
199  G4String titl = title[ih];
200  G4double vmin = valmin, vmax = valmax;
201  Unit[ih] = 1.;
202 
203  if (unit != "none") {
204  titl = title[ih] + " (" + unit + ")";
205  Unit[ih] = G4UnitDefinition::GetValueOf(unit);
206  vmin = valmin/Unit[ih]; vmax = valmax/Unit[ih];
207  }
208 
209  exist[ih] = true;
210  Label[ih] = id[ih];
211  Title[ih] = titl;
212  Nbins[ih] = nbins;
213  Vmin[ih] = vmin;
214  Vmax[ih] = vmax;
215  Width[ih] = (valmax-valmin)/nbins;
216 
217  G4cout << "----> SetHisto " << ih << ": " << titl << "; "
218  << nbins << " bins from "
219  << vmin << " " << unit << " to " << vmax << " " << unit << G4endl;
220 
221 }
222 
223 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
224 
226 {
227  if (ih >= MaxHisto) {
228  G4cout << "---> warning from HistoManager::Normalize() : histo " << ih
229  << " fac= " << fac << G4endl;
230  return;
231  }
232 #ifdef G4ANALYSIS_USE
233  if(exist[ih]) histo[ih]->scale(fac);
234 #endif
235 }
236 
237 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
238 
240 {
241  if (ih >= MaxHisto) {
242  G4cout << "---> warning from HistoManager::RemoveHisto() : histo " << ih
243  << "does not exist" << G4endl;
244  return;
245  }
246 
247  histo[ih] = 0; exist[ih] = false;
248 }
249 
250 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
251 
253 {
254  if (ih < MaxHisto) { ascii[ih] = true; ascii[0] = true; }
255  else
256  G4cout << "---> warning from HistoManager::PrintHisto() : histo " << ih
257  << "does not exist" << G4endl;
258 }
259 
260 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
261 
262 #include <fstream>
263 
264 void HistoManager::saveAscii()
265 {
266 #ifdef G4ANALYSIS_USE
267 
268  if (!ascii[0]) return;
269 
270  G4String name = fileName[0] + ".ascii";
271  std::ofstream File(name, std::ios::out);
272  File.setf( std::ios::scientific, std::ios::floatfield );
273 
274  //write selected histograms
275  for (G4int ih=0; ih<MaxHisto; ih++) {
276  if (exist[ih] && ascii[ih]) {
277  File << "\n 1D histogram " << ih << ": " << Title[ih]
278  << "\n \n \t X \t\t Y" << G4endl;
279 
280  for (G4int iBin=0; iBin<Nbins[ih]; iBin++) {
281  File << " " << iBin << "\t"
282  << 0.5*(histo[ih]->axis().binLowerEdge(iBin) +
283  histo[ih]->axis().binUpperEdge(iBin)) << "\t"
284  << histo[ih]->binHeight(iBin)
285  << G4endl;
286  }
287  }
288  }
289 #endif
290 }
291 
292 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
293 
295 {
296  if (nt < 1 || nt >= MaxNtupl) {
297  G4cout << "---> warning from HistoManager::SetNtuple() : Ntuple " << nt
298  << "does not exist" << G4endl;
299  return;
300  }
301 
302  const G4String id[] = { "100", "101" };
303  const G4String title[] =
304  { "dummy", //0
305  "Energy deposit in subModule" //1
306  };
307 
308  G4String column[MaxNtupl];
309 
310  column[0] = " int dum=0 "; //0
311 
312  column[1] =
313  "double Evis0, Evis1, Evis2, Evis3, Evis4, Evis5, Evis6, Evis7, Evis8, Evis9";
314  column[1] +=
315  ", Evis10, Evis11, Evis12, Evis13, Evis14, Evis15, Evis16, Evis17";
316  column[1] +=
317  ", Edep0, Edep1, Edep2, Edep3, Edep4, Edep5, Edep6, Edep7, Edep8, Edep9";
318  column[1] +=
319  ", Edep10, Edep11, Edep12, Edep13, Edep14, Edep15, Edep16, Edep17";
320 
321  G4String titl = title[nt];
322 
323  existNt[nt] = true;
324  LabelNt[nt] = id[nt];
325  TitleNt[nt] = titl;
326  ColumnNt[nt] = column[nt];
327 
328  G4cout << "----> SetNtuple " << nt << ": " << titl << "; " << G4endl;
329 
330 }
331 
332 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
333 
335 {
336  if (nt >= MaxNtupl) {
337  G4cout << "---> warning from HistoManager::FillNtuple() : Ntuple " << nt
338  << " does not exist " << column << value << G4endl;
339  return;
340  }
341 #ifdef G4ANALYSIS_USE
342  if(existNt[nt]) ntupl[nt]->fill(column, value);
343 #endif
344 }
345 
346 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
347 
349 {
350  if (nt >= MaxNtupl) {
351  G4cout << "---> warning from HistoManager::AddRowNtuple() : Ntuple " << nt
352  << " do not exist" << G4endl;
353  return;
354  }
355 #ifdef G4ANALYSIS_USE
356  if(existNt[nt]) ntupl[nt]->addRow();
357 #endif
358 }
359 
360 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
361 
362 
void FillNtuple(G4int id, G4int column, G4double value)
void AddRowNtuple(G4int id)
void FillHisto(G4int id, G4double bin, G4double weight=1.0)
const XML_Char * name
Definition: expat.h:151
double weight
Definition: plottest35.C:25
int G4int
Definition: G4Types.hh:78
void Normalize(G4int id, G4double fac)
static G4double GetValueOf(const G4String &)
G4GLOB_DLL std::ostream G4cout
tuple tree
Definition: gammaraytel.py:4
bool G4bool
Definition: G4Types.hh:79
TTree * nt
Definition: plotHisto.C:21
const G4int MaxNtupl
Definition: HistoManager.hh:47
void SetHisto(G4int, G4int, G4double, G4double, const G4String &unit="none")
subroutine title(NA, NB, NCA, NCB)
Definition: dpm25nuc7.f:1744
const G4int MaxHisto
Definition: HistoManager.hh:47
const XML_Char int const XML_Char * value
Definition: expat.h:331
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
void RemoveHisto(G4int)
void PrintHisto(G4int)
void SetNtuple(G4int)
tuple hf
Definition: gammaraytel.py:8