source: trunk/examples/extended/radioactivedecay/exrdm/src/exrdmHisto.cc @ 850

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

update
File size: 11.4 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#include "exrdmHisto.hh"
28#include "exrdmHistoMessenger.hh"
29#include "G4ParticleTable.hh"
30
31#include "G4Tokenizer.hh"
32
33#ifdef G4ANALYSIS_USE_AIDA
34#include <AIDA/AIDA.h>
35#endif
36//
37#ifdef G4ANALYSIS_USE_ROOT
38#include "TROOT.h"
39#include "TApplication.h"
40#include "TGClient.h"
41#include "TCanvas.h"
42#include "TSystem.h"
43#include "TTree.h"
44#include "TBranch.h"
45#include "TFile.h"
46#include "TH1D.h"
47#include "TNtuple.h"
48#endif
49//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
50exrdmHisto::exrdmHisto()
51{
52  verbose    = 0;
53  histName   = "exrdm";
54  histType   = "root";
55  nHisto     = 0;
56  nTuple     = 0;
57  defaultAct = 1;
58  //
59#ifdef G4ANALYSIS_USE_AIDA
60  histo.clear();
61  ntup.clear();
62#endif
63#ifdef G4ANALYSIS_USE_ROOT
64  ROOThisto.clear();
65  ROOTntup.clear();
66  Rarray.clear();
67  Rcol.clear();
68#endif
69  active.clear();
70  bins.clear();
71  xmin.clear();
72  xmax.clear();
73  unit.clear();
74  ids.clear();
75  titles.clear();
76  tupleName.clear();
77  tupleId.clear();
78  tupleList.clear();
79  tupleListROOT.clear();
80  messenger  = 0;
81
82  messenger = new exrdmHistoMessenger(this);
83}
84
85//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
86
87exrdmHisto::~exrdmHisto()
88{
89#ifdef G4ANALYSIS_USE_AIDA
90  for(G4int i=0; i<nHisto; i++) {
91    if(histo[i]) delete histo[i];
92  }
93#endif
94#ifdef G4ANALYSIS_USE_ROOT
95  for(G4int i=0; i<nHisto; i++) {
96    if(ROOThisto[i]) delete ROOThisto[i];
97  }
98#endif
99  delete messenger;
100}
101
102//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
103
104void exrdmHisto::book()
105{
106  G4cout << "### exrdmHisto books " << nHisto << " histograms " << G4endl; 
107#ifdef G4ANALYSIS_USE_AIDA
108  // Creating the analysis factory
109  AIDA::IAnalysisFactory* af = AIDA_createAnalysisFactory();
110  // Creating the tree factory
111  AIDA::ITreeFactory* tf = af->createTreeFactory(); 
112  // Creating a tree mapped to a new aida file.
113  G4String aidaFileName;
114  if (fileType == "hbook")
115    aidaFileName = histName +G4String(".hbook");
116  else
117    aidaFileName = histName +G4String(".aida");
118  tree = tf->create(aidaFileName,histType,false,true,"uncompress");
119  if(tree) 
120    G4cout << "Tree store  : " << tree->storeName() << G4endl;
121  else
122    G4cout << "ERROR: Tree store " << histName  << " is not created!" << G4endl;
123
124  // Creating a histogram factory, whose histograms will be handled by the tree
125  AIDA::IHistogramFactory* hf = af->createHistogramFactory(*tree);
126  // Creating an 1-dimensional histograms in the root directory of the tree
127  for(G4int i=0; i<nHisto; i++) {
128    if(active[i]) {
129      histo[i] = hf->createHistogram1D(ids[i], titles[i], bins[i], xmin[i], xmax[i]);
130    }
131  }
132  // Creating a tuple factory, whose tuples will be handled by the tree 
133  AIDA::ITupleFactory* tpf =  af->createTupleFactory( *tree );
134  for(G4int i=0; i<nTuple; i++) {
135    if(tupleList[i] != "") {
136      G4cout << "Creating Ntuple: " << tupleName[i] << G4endl;
137      ntup[i] = tpf->create(tupleId[i], tupleName[i], tupleList[i]);
138    }
139  }
140#endif
141
142#ifdef G4ANALYSIS_USE_ROOT
143  new TApplication("App", ((int *)0), ((char **)0));
144  G4String fileNameROOT = histName + G4String(".root");
145  hfileROOT = new TFile(fileNameROOT.c_str() ,"RECREATE","ROOT file for exRDM");
146 
147  // Creating an 1-dimensional histograms in the root directory of the tree
148  for(G4int i=0; i<nHisto; i++) {
149    if(active[i]) {
150      G4String id = G4String("h")+ids[i];
151      ROOThisto[i] = new TH1D(id, titles[i], bins[i], xmin[i], xmax[i]);
152      G4cout << "ROOT Histo " << ids[i] << " " << titles[i] << " booked " << G4endl;
153    }
154  }
155  // Now the ntuples 
156  for(G4int i=0; i<nTuple; i++) {
157    if(tupleListROOT[i] != "") {
158      G4String id = G4String("t")+tupleId[i];
159      G4cout << "Creating Ntuple "<<tupleId[i] << " in ROOT file: " << tupleName[i] << G4endl;
160      ROOTntup[i] = new TNtuple(id, tupleName[i], tupleListROOT[i]);
161    }
162  }
163#endif
164
165} 
166
167//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
168
169void exrdmHisto::save()
170{
171#ifdef G4ANALYSIS_USE_AIDA
172  // Write histogram file
173  tree->commit();
174  G4cout << "Closing the AIDA tree..." << G4endl;
175  tree->close();
176  G4cout << "Histograms and Ntuples are saved" << G4endl;
177#endif
178#ifdef G4ANALYSIS_USE_ROOT
179  G4cout << "ROOT: files writing..." << G4endl;
180  hfileROOT->Write();
181  G4cout << "ROOT: files closing..." << G4endl;
182  hfileROOT->Close();
183#endif
184}
185
186
187//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
188
189void exrdmHisto::add1D(const G4String& id, const G4String& name, G4int nb, 
190                  G4double x1, G4double x2, G4double u)
191{
192  if(verbose > 0) {
193    G4cout << "New histogram will be booked: #" << id << "  <" << name
194           << "  " << nb << "  " << x1 << "  " << x2 << "  " << u
195           << G4endl;
196  }
197  nHisto++;
198  x1 /= u;
199  x2 /= u;
200  active.push_back(defaultAct);
201  bins.push_back(nb);
202  xmin.push_back(x1);
203  xmax.push_back(x2);
204  unit.push_back(u);
205  ids.push_back(id);
206  titles.push_back(name);
207#ifdef G4ANALYSIS_USE_AIDA
208  histo.push_back(0);
209#endif
210#ifdef G4ANALYSIS_USE_ROOT
211  ROOThisto.push_back(0);
212#endif
213}
214
215//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
216
217void exrdmHisto::setHisto1D(G4int i, G4int nb, G4double x1, G4double x2, G4double u)
218{
219  if(i>=0 && i<nHisto) {
220    if(verbose > 0) {
221      G4cout << "Update histogram: #" << i 
222             << "  " << nb << "  " << x1 << "  " << x2 << "  " << u
223             << G4endl;
224    }
225    bins[i] = nb;
226    xmin[i] = x1;
227    xmax[i] = x2;
228    unit[i] = u;
229  } else {
230    G4cout << "exrdmHisto::setexrdmHisto1D: WARNING! wrong histogram index " << i << G4endl;
231  }
232}
233
234//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
235
236void exrdmHisto::fillHisto(G4int i, G4double x, G4double w)
237{
238  if(verbose > 1) {
239    G4cout << "fill histogram: #" << i << " at x= " << x
240           << "  weight= " << w
241           << G4endl;   
242  }
243#ifdef G4ANALYSIS_USE_AIDA 
244  if(i>=0 && i<nHisto) {
245    histo[i]->fill((float)(x/unit[i]), (float)w);
246  } else {
247    G4cout << "exrdmHisto::fill: WARNING! wrong AIDA histogram index " << i << G4endl;
248  }
249#endif
250#ifdef G4ANALYSIS_USE_ROOT 
251  if(i>=0 && i<nHisto) {
252    ROOThisto[i]->Fill(x/unit[i],w);
253  } else {
254    G4cout << "exrdmHisto::fill: WARNING! wrong ROOT histogram index " << i << G4endl;
255  }
256#endif
257}
258
259//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
260
261void exrdmHisto::scaleHisto(G4int i, G4double x)
262{
263  if(verbose > 0) {
264    G4cout << "Scale histogram: #" << i << " by factor " << x << G4endl;   
265  }
266#ifdef G4ANALYSIS_USE_AIDA 
267  if(i>=0 && i<nHisto) {
268    histo[i]->scale(x);
269    G4cout << "exrdmHisto::scale: WARNING! wrong AIDA histogram index " << i << G4endl;
270  }
271#endif
272#ifdef G4ANALYSIS_USE_ROOT 
273  if(i>=0 && i<nHisto) {
274    ROOThisto[i]->Scale(x);
275  } else {
276    G4cout << "exrdmHisto::scale: WARNING! wrong ROOT histogram index " << i << G4endl;
277  }
278#endif
279}
280
281//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
282
283void exrdmHisto::addTuple(const G4String& w1, const G4String& w2, const G4String& w3)
284{
285//  G4cout << w1 << " " << w2 << " " << w3 << G4endl;
286  std::vector<float> ar;
287  ar.clear();
288  for (size_t i = 0; i < 20; i++) ar.push_back(0.);
289  nTuple++;
290#ifdef G4ANALYSIS_USE_ROOT
291  Rarray.push_back(ar);
292#endif
293  tupleId.push_back(w1);
294  tupleName.push_back(w2) ;
295  // convert AIDA header to ROOT header for ntuple
296  G4Tokenizer next(w3);
297  G4String token = next();
298  G4String ROOTList1 = "" ;
299  G4int col = 0;
300  while ( token != "") {
301   token = next();
302   ROOTList1 = ROOTList1 + token +G4String(":");
303   col++;
304  }
305  G4String ROOTList = ROOTList1.substr(0,ROOTList1.length()-2);
306//  G4cout << ROOTList << G4endl;
307  tupleListROOT.push_back(ROOTList);
308  //
309#ifdef G4ANALYSIS_USE_AIDA 
310  ntup.push_back(0);
311#endif
312#ifdef G4ANALYSIS_USE_ROOT
313  ROOTntup.push_back(0);
314  Rcol.push_back(col-1);
315#endif
316}
317
318//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
319
320void exrdmHisto::fillTuple(G4int i, const G4String& parname, G4double x)
321{
322  if(verbose > 1) 
323    G4cout << "fill tuple # " << i
324           <<" with  parameter <" << parname << "> = " << x << G4endl; 
325#ifdef G4ANALYSIS_USE_AIDA 
326  if(ntup[i]) ntup[i]->fill(ntup[i]->findColumn(parname), (float)x);
327#endif
328}
329
330//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
331
332void exrdmHisto::fillTuple(G4int i, G4int col, G4double x)
333{
334  if(verbose > 1) {
335    G4cout << "fill tuple # " << i
336           <<" in column < " << col << "> = " << x << G4endl; 
337  }
338#ifdef G4ANALYSIS_USE_AIDA 
339  if(ntup[i]) ntup[i]->fill(col, (float)x);
340#endif
341
342#ifdef G4ANALYSIS_USE_ROOT 
343  if(ROOTntup[i]) (Rarray[i])[col] = float(x);
344#endif
345
346}
347
348//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
349
350void exrdmHisto::fillTuple(G4int i, const G4String& parname, G4String x)
351{
352  if(verbose > 1) {
353    G4cout << "fill tuple # " << i
354           <<" with  parameter <" << parname << "> = " << x << G4endl; 
355  }
356#ifdef G4ANALYSIS_USE_AIDA 
357  if(ntup[i]) ntup[i]->fill(ntup[i]->findColumn(parname), x);
358#endif
359}
360
361//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
362
363void exrdmHisto::addRow(G4int i)
364{
365#ifdef G4ANALYSIS_USE_AIDA
366  if(ntup[i]) ntup[i]->addRow();
367#endif
368
369#ifdef G4ANALYSIS_USE_ROOT
370  float ar[4];
371  for (G4int j=0; j < Rcol[i]; j++) {
372//      G4cout << i << " " << Rarray[i][j] << G4endl;
373      ar[j] = Rarray[i][j];       
374  } 
375  if(ROOTntup[i]) ROOTntup[i]->Fill(ar);
376#endif
377
378} 
379
380//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
381
382void exrdmHisto::setFileName(const G4String& nam) 
383{
384  histName = nam;
385}
386
387//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
388
389void exrdmHisto::setFileType(const G4String& nam) 
390{
391  histType = nam;
392}
393
394//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
395
396const G4String& exrdmHisto::FileType() const
397{
398  return histType;
399}
400
401//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
402
Note: See TracBrowser for help on using the repository browser.