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 //
28 //
29 // $Id: HistoManager.cc 73005 2013-08-15 08:14:21Z gcosmo $
30 //
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
33 
34 #include "HistoManager.hh"
35 #include "HistoMessenger.hh"
36 #include "G4UnitsTable.hh"
37 
38 #ifdef G4ANALYSIS_USE
39 #include "AIDA/AIDA.h"
40 #endif
41 
42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
43 
45 #ifdef G4ANALYSIS_USE
46 :af(0),tree(0),factoryOn(false)
47 #endif
48 {
49 #ifdef G4ANALYSIS_USE
50  // Creating the analysis factory
51  af = AIDA_createAnalysisFactory();
52  if(!af) {
53  G4cout << " HistoManager::HistoManager() :"
54  << " problem creating the AIDA analysis factory."
55  << G4endl;
56  }
57 #endif
58 
59  fileName[0] = "Pol01.aida";
60  fileType = "xml";
61  fileOption = "";
62  // histograms
63  for (G4int k=0; k<MaxHisto; k++) {
64  histo[k] = 0;
65  exist[k] = false;
66  Unit[k] = 1.0;
67  Width[k] = 1.0;
68  }
69 
70  histoMessenger = new HistoMessenger(this);
71 }
72 
73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
74 
76 {
77  delete histoMessenger;
78 
79 #ifdef G4ANALYSIS_USE
80  delete af;
81 #endif
82 }
83 
84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
85 
86 void HistoManager::book()
87 {
88 #ifdef G4ANALYSIS_USE
89  if(!af) return;
90 
91  // Creating a tree mapped to an aida file.
92  if (fileName[0].find(".")==G4String::npos) {
93  fileName[1] = fileName[0] + "." + fileType;
94  }
95  else {
96  fileName[1] = fileName[0];
97  }
98 
99  G4bool readOnly = false;
100  G4bool createNew = true;
101  AIDA::ITreeFactory* tf = af->createTreeFactory();
102  tree = tf->create(fileName[1], fileType, readOnly, createNew, fileOption);
103  delete tf;
104  if(!tree) {
105  G4cout << "HistoManager::book() :"
106  << " problem creating the AIDA tree with "
107  << " storeName = " << fileName[1]
108  << " storeType = " << fileType
109  << " readOnly = " << readOnly
110  << " createNew = " << createNew
111  << " options = " << fileOption
112  << G4endl;
113  return;
114  }
115 
116  // Creating a histogram factory, whose histograms will be handled by the tree
117  AIDA::IHistogramFactory* hf = af->createHistogramFactory(*tree);
118 
119  // create selected histograms
120  for (G4int k=0; k<MaxHisto; k++) {
121  if (exist[k]) {
122  histo[k] = hf->createHistogram1D( Label[k], Title[k],
123  Nbins[k], Vmin[k], Vmax[k]);
124  factoryOn = true;
125  }
126  }
127  delete hf;
128  if(factoryOn)
129  G4cout << "\n----> Histogram Tree is opened " << fileName[1] << G4endl;
130 #endif
131 }
132 
133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
134 
135 void HistoManager::save()
136 {
137 #ifdef G4ANALYSIS_USE
138  if (factoryOn) {
139  tree->commit(); // Writing the histograms to the file
140  tree->close(); // and closing the tree (and the file)
141  G4cout << "\n----> Histogram Tree is saved in " << fileName[1] << G4endl;
142 
143  delete tree;
144  tree = 0;
145  factoryOn = false;
146  }
147 #endif
148 }
149 
150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
151 
152 void HistoManager::FillHistos(const G4String & particleName,
153  G4double kinEnergy, G4double costheta,
154  G4double phi,
155  G4double longitudinalPolarization)
156 {
157  G4int id=1;
158  if (particleName=="gamma") id=1;
159  else if (particleName=="e-") id=5;
160  else if (particleName=="e+") id=9;
161  else return;
162 
163  if (costheta>=1.) costheta=.99999999;
164  if (costheta<-1.) costheta=-1.;
165  FillHisto(id,kinEnergy);
166  FillHisto(id+1,costheta);
167  FillHisto(id+2,phi);
168  FillHisto(id+3,longitudinalPolarization);
169 }
170 
172 {
173  if (ih > MaxHisto) {
174  G4cout << "---> warning from HistoManager::FillHisto() : histo " << ih
175  << "does not exist; e= " << e << " w= " << weight << G4endl;
176  return;
177  }
178 #ifdef G4ANALYSIS_USE
179  if(exist[ih]) histo[ih]->fill(e/Unit[ih], weight);
180 #endif
181 }
182 
183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
184 
185 void HistoManager::SetHisto(G4int ih, G4int nbins, G4double valmin,
186  G4double valmax, const G4String& unit)
187 {
188  if (ih > MaxHisto) {
189  G4cout << "---> warning from HistoManager::SetHisto() : histo " << ih
190  << "does not exist" << G4endl;
191  return;
192  }
193 
194  const G4String id[] = { "0", "1", "2", "3", "4", "5",
195  "6", "7", "8", "9", "10", "11", "12"};
196  const G4String title[] =
197  { "dummy", //0
198  "Gamma Energy distribution", //1
199  "Gamma Cos(Theta) distribution", //2
200  "Gamma Phi angular distribution", //3
201  "Gamma longitudinal Polarization", //4
202  "Electron Energy distribution", //5
203  "Electron Cos(Theta) distribution", //6
204  "Electron Phi angular distribution", //7
205  "Electron longitudinal Polarization", //8
206  "Positron Energy distribution", //9
207  "Positron Cos(Theta) distribution", //10
208  "Positron Phi angular distribution", //11
209  "Positron longitudinal Polarization" //12
210  };
211 
212 
213  G4String titl = title[ih];
214  G4double vmin = valmin, vmax = valmax;
215  Unit[ih] = 1.;
216 
217  if (unit != "none") {
218  titl = title[ih] + " (" + unit + ")";
219  Unit[ih] = G4UnitDefinition::GetValueOf(unit);
220  vmin = valmin/Unit[ih]; vmax = valmax/Unit[ih];
221  }
222 
223  exist[ih] = true;
224  Label[ih] = id[ih];
225  Title[ih] = titl;
226  Nbins[ih] = nbins;
227  Vmin[ih] = vmin;
228  Vmax[ih] = vmax;
229  Width[ih] = (valmax-valmin)/nbins;
230 
231  G4cout << "----> SetHisto " << ih << ": " << titl << "; "
232  << nbins << " bins from "
233  << vmin << " " << unit << " to " << vmax << " " << unit << G4endl;
234 
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 
void FillHisto(G4int id, G4double bin, G4double weight=1.0)
double weight
Definition: plottest35.C:25
int G4int
Definition: G4Types.hh:78
void FillHistos(const G4String &particleName, G4double kinEnergy, G4double costheta, G4double phi, G4double longitudinalPolarization)
static G4double GetValueOf(const G4String &)
G4GLOB_DLL std::ostream G4cout
tuple tree
Definition: gammaraytel.py:4
bool G4bool
Definition: G4Types.hh:79
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
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
void RemoveHisto(G4int)
tuple hf
Definition: gammaraytel.py:8