source: trunk/examples/extended/polarisation/Pol01/src/HistoManager.cc @ 850

Last change on this file since 850 was 807, checked in by garnier, 16 years ago

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