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

Last change on this file since 1036 was 807, checked in by garnier, 17 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.