Geant4  10.01.p01
exrdmHisto.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: exrdmHisto.cc 78592 2014-01-08 10:30:37Z gcosmo $
27 //
30 //
31 #ifdef G4ANALYSIS_USE
32 #include <AIDA/AIDA.h>
33 #endif
34 //
35 #ifdef G4ANALYSIS_USE_ROOT
36 #include "TROOT.h"
37 #include "TApplication.h"
38 #include "TGClient.h"
39 #include "TCanvas.h"
40 #include "TSystem.h"
41 #include "TTree.h"
42 #include "TBranch.h"
43 #include "TFile.h"
44 #include "TH1D.h"
45 #include "TNtuple.h"
46 #endif
47 
48 #include "exrdmHisto.hh"
49 #include "exrdmHistoMessenger.hh"
50 #include "G4ParticleTable.hh"
51 
52 #include "G4Tokenizer.hh"
53 
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
56 :fHistName("exrdm"), fHistType("root"),
57  fNHisto(0), fNTuple(0), fVerbose(0),
58  fDefaultAct(1)
59 {
60 #ifdef G4ANALYSIS_USE
61  fAida = 0;
62  fTree = 0;
63 #endif
64 
65 #ifdef G4ANALYSIS_USE_ROOT
66  fROOThisto.clear();
67  fROOTntup.clear();
68  fRarray.clear();
69  fRcol.clear();
70 #endif
71 
72  fActive.clear();
73  fBins.clear();
74  fXmin.clear();
75  fXmax.clear();
76  fUnit.clear();
77  fIds.clear();
78  fTitles.clear();
79  fTupleName.clear();
80  fTupleId.clear();
81  fTupleList.clear();
82  fTupleListROOT.clear();
83 
84  fMessenger = new exrdmHistoMessenger(this);
85 }
86 
87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
88 
90 {
91 #ifdef G4ANALYSIS_USE
92  fHisto.clear();
93  fNtup.clear();
94 #endif
95 #ifdef G4ANALYSIS_USE_ROOT
96  //FIXME : G.Barrand : the below is crashy.
97  // In principle the TH are deleted
98  // when doing the TFile::Close !
99  // In fact the fHfileROOT should
100  // be deleted in Save(). And I am pretty
101  // sure that the TApplication is not needed.
102  //
103  // removed by F.Lei
104  // for(G4int i=0; i<fNHisto; i++) {
105  // if(fROOThisto[i]) delete fROOThisto[i];
106  // }
107 #endif
108  delete fMessenger;
109 }
110 
111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
112 
114 {
115 #ifdef G4ANALYSIS_USE
116  G4cout << "### exrdmHisto books " << fNHisto << " histograms " << G4endl;
117  // Creating the analysis factory
118  fAida = AIDA_createAnalysisFactory();
119  if(!fAida) {
120  G4cout << "ERROR: can't get AIDA." << G4endl;
121  return;
122  }
123  // Creating the fTree factory
124  {AIDA::ITreeFactory* tf = fAida->createTreeFactory();
125  // Creating a fTree mapped to a new fAida file.
126  G4String fileName = fHistName + "." + fHistType;
127  if (fHistType == "root") fileName = fHistName + "_aida." + fHistType;
128  fTree = tf->create(fileName,fHistType,false,true,"");
129  delete tf;
130  if(!fTree) {
131  G4cout << "ERROR: Tree store " << fHistName << " is not created!" << G4endl;
132  return;
133  }
134  G4cout << "Tree store : " << fTree->storeName() << G4endl;}
135  // Creating a histogram factory, whose histograms will be handled by the fTree
136  {AIDA::IHistogramFactory* hf = fAida->createHistogramFactory(*fTree);
137  // Creating an 1-dimensional histograms in the root directory of the fTree
138  for(G4int i=0; i<fNHisto; i++) {
139  if(fActive[i]) {
140  if(fVerbose>1)
141  G4cout<<"Book: histogram "<< i << " id= " << fIds[i] <<G4endl;
142  G4String tit = fIds[i];
143  if(fHistType == "root") tit = "h" + fIds[i];
144  fHisto[i] = hf->createHistogram1D(tit, fTitles[i], fBins[i], fXmin[i],
145  fXmax[i]);
146  }
147  }
148  delete hf;
149  G4cout << "AIDA histograms are booked" << G4endl;}
150 
151  // Creating a tuple factory, whose tuples will be handled by the fTree
152  {AIDA::ITupleFactory* tpf = fAida->createTupleFactory( *fTree );
153  G4cout << "AIDA will Book " << fNTuple << " ntuples" << G4endl;
154  for(G4int i=0; i<fNTuple; i++) {
155  if(fTupleList[i] != "") {
156  G4cout << "Creating Ntuple: " << fTupleName[i] <<":" <<fTupleList[i]
157  << G4endl;
158  fNtup[i] = tpf->create(fTupleId[i], fTupleName[i], fTupleList[i],"");
159  }
160  }
161  delete tpf;
162  G4cout << "AIDA ntuples are booked" << G4endl;}
163 #endif
164 
165 #ifdef G4ANALYSIS_USE_ROOT
166 // new TApplication("App", ((int *)0), ((char **)0));
167  G4String fileNameROOT = fHistName + G4String(".root");
168  fHfileROOT = new TFile(fileNameROOT.c_str() ,"RECREATE","ROOT file for exRDM");
169  G4cout << "Root file: " << fileNameROOT << G4endl;
170  // Creating an 1-dimensional histograms in the root directory of the fTree
171  for(G4int i=0; i<fNHisto; i++) {
172  if(fActive[i]) {
173  G4String id = G4String("h")+fIds[i];
174  fROOThisto[i] = new TH1D(id, fTitles[i], fBins[i], fXmin[i], fXmax[i]);
175  G4cout << "ROOT Histo " << fIds[i] << " " << fTitles[i] << " booked "
176  << G4endl;
177  }
178  }
179  // Now the ntuples
180  for(G4int i=0; i<fNTuple; i++) {
181  if(fTupleListROOT[i] != "") {
182  G4String id = G4String("t")+fTupleId[i];
183  G4cout << "Creating Ntuple "<<fTupleId[i] << " in ROOT file: "
184  << fTupleName[i] << G4endl;
185  fROOTntup[i] = new TNtuple(id, fTupleName[i], fTupleListROOT[i]);
186  G4cout << "ROOT Ntuple " << id << " " << fTupleName[i] <<" "
187  << fTupleListROOT[i]<< " booked " << G4endl;
188  }
189  }
190 #endif
191 
192 }
193 
194 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
195 
197 {
198 #ifdef G4ANALYSIS_USE
199  // Write histogram file
200  fTree->commit();
201  G4cout << "Closing the AIDA fTree..." << G4endl;
202  fTree->close();
203  G4cout << "Histograms and Ntuples are saved" << G4endl;
204  delete fTree;
205  fTree = 0;
206  delete fAida;
207  fAida = 0;
208  {for(G4int i=0; i<fNHisto; i++) fHisto[i] = 0;}
209  {for(G4int i=0; i<fNTuple; i++) fNtup[i] = 0;}
210 #endif
211 #ifdef G4ANALYSIS_USE_ROOT
212  G4cout << "ROOT: files writing..." << G4endl;
213  fHfileROOT->Write();
214  G4cout << "ROOT: files closing..." << G4endl;
215  fHfileROOT->Close();
216  //
217  // F.Lei added following Guy's suggestion!
218  delete fHfileROOT;
219 
220 #endif
221 }
222 
223 
224 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
225 
226 void exrdmHisto::Add1D(const G4String& id, const G4String& name, G4int nb,
227  G4double x1, G4double x2, G4double u)
228 {
229  if(fVerbose > 0) {
230  G4cout << "New histogram will be booked: #" << id << " <" << name
231  << " " << nb << " " << x1 << " " << x2 << " " << u
232  << G4endl;
233  }
234  fNHisto++;
235  x1 /= u;
236  x2 /= u;
237  fActive.push_back(fDefaultAct);
238  fBins.push_back(nb);
239  fXmin.push_back(x1);
240  fXmax.push_back(x2);
241  fUnit.push_back(u);
242  fIds.push_back(id);
243  fTitles.push_back(name);
244 #ifdef G4ANALYSIS_USE
245  fHisto.push_back(0);
246 #endif
247 #ifdef G4ANALYSIS_USE_ROOT
248  fROOThisto.push_back(0);
249 #endif
250 }
251 
252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
253 
255 {
256  if(i>=0 && i<fNHisto) {
257  if(fVerbose > 0) {
258  G4cout << "Update histogram: #" << i
259  << " " << nb << " " << x1 << " " << x2 << " " << u
260  << G4endl;
261  }
262  fBins[i] = nb;
263  fXmin[i] = x1;
264  fXmax[i] = x2;
265  fUnit[i] = u;
266  } else {
267  G4cout << "exrdmHisto::setexrdmHisto1D: WARNING! wrong histogram index "
268  << i << G4endl;
269  }
270 }
271 
272 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
273 
275 {
276  if(fVerbose > 1) {
277  G4cout << "fill histogram: #" << i << " at x= " << x
278  << " weight= " << w
279  << G4endl;
280  }
281 #ifdef G4ANALYSIS_USE
282  if(i>=0 && i<fNHisto) {
283  fHisto[i]->fill(x/fUnit[i], w);
284  } else {
285  G4cout << "exrdmHisto::fill: WARNING! wrong AIDA histogram index "
286  << i << G4endl;
287  }
288 #endif
289 #ifdef G4ANALYSIS_USE_ROOT
290  if(i>=0 && i<fNHisto) {
291  fROOThisto[i]->Fill(x/fUnit[i],w);
292  } else {
293  G4cout << "exrdmHisto::fill: WARNING! wrong ROOT histogram index "
294  << i << G4endl;
295  }
296 #endif
297 }
298 
299 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
300 
302 {
303  if(fVerbose > 0) {
304  G4cout << "Scale histogram: #" << i << " by factor " << x << G4endl;
305  }
306 #ifdef G4ANALYSIS_USE
307  if(i>=0 && i<fNHisto) {
308  fHisto[i]->scale(x);
309  G4cout << "exrdmHisto::scale: WARNING! wrong AIDA histogram index "
310  << i << G4endl;
311  }
312 #endif
313 #ifdef G4ANALYSIS_USE_ROOT
314  if(i>=0 && i<fNHisto) {
315  fROOThisto[i]->Scale(x);
316  } else {
317  G4cout << "exrdmHisto::scale: WARNING! wrong ROOT histogram index "
318  << i << G4endl;
319  }
320 #endif
321 }
322 
323 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
324 
325 #ifdef G4ANALYSIS_USE
326 void exrdmHisto::AddTuple(const G4String& w1, const G4String& w2,
327  const G4String& w3 )
328 #else
329 #ifdef G4ANALYSIS_USE_ROOT
330 void exrdmHisto::AddTuple(const G4String& w1, const G4String& w2,
331  const G4String& w3 )
332 #else
333 void exrdmHisto::AddTuple(const G4String& w1, const G4String& w2,
334  const G4String& )
335 #endif
336 #endif
337 
338 {
339  //G4cout << w1 << " " << w2 << " " << w3 << G4endl;
340  fNTuple++;
341  fTupleId.push_back(w1);
342  fTupleName.push_back(w2) ;
343 #ifdef G4ANALYSIS_USE
344  fTupleList.push_back(w3);
345  fNtup.push_back(0);
346 #endif
347 
348 #ifdef G4ANALYSIS_USE_ROOT
349  std::vector<float> ar;
350  ar.clear();
351  for (size_t i = 0; i < 20; i++) ar.push_back(0.);
352  fRarray.push_back(ar);
353  // convert AIDA header to ROOT header for ntuple
354  G4Tokenizer next(w3);
355  G4String token = next();
356  G4String ROOTList1 = "" ;
357  G4int col = 0;
358  while ( token != "") {
359  token = next();
360  if (token == ",") token = next();
361  if (token.contains(",")) token.remove(token.size()-1);
362  ROOTList1 = ROOTList1 + token + G4String(":");
363  col++;
364  }
365  G4String ROOTList = ROOTList1.substr(0,ROOTList1.length()-2);
366 // G4cout << ROOTList << G4endl;
367  fTupleListROOT.push_back(ROOTList);
368  fROOTntup.push_back(0);
369  fRcol.push_back(col-1);
370 #endif
371 }
372 
373 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
374 
375 void exrdmHisto::FillTuple(G4int i, const G4String& parname, G4double x)
376 {
377  if(fVerbose > 1)
378  G4cout << "fill tuple # " << i
379  <<" with parameter <" << parname << "> = " << x << G4endl;
380 #ifdef G4ANALYSIS_USE
381  if(fNtup[i]) fNtup[i]->fill(fNtup[i]->findColumn(parname), x);
382 #endif
383 }
384 
385 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
386 
388 {
389  if(fVerbose > 1) {
390  G4cout << "fill tuple # " << i
391  <<" in column < " << col << "> = " << x << G4endl;
392  }
393 #ifdef G4ANALYSIS_USE
394  if(fNtup[i]) fNtup[i]->fill(col,double(x));
395 #endif
396 
397 #ifdef G4ANALYSIS_USE_ROOT
398  if(fROOTntup[i]) (fRarray[i])[col] = float(x);
399 #endif
400 
401 }
402 
403 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
404 
405 void exrdmHisto::FillTuple(G4int i, const G4String& parname, G4String& x)
406 {
407  if(fVerbose > 1) {
408  G4cout << "fill tuple # " << i
409  <<" with parameter <" << parname << "> = " << x << G4endl;
410  }
411 #ifdef G4ANALYSIS_USE
412  if(fNtup[i]) fNtup[i]->fill(fNtup[i]->findColumn(parname), x);
413 #endif
414 
415 }
416 
417 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
418 
420 {
421  if(fVerbose > 1) G4cout << "Added a raw #" << i << " to tuple" << G4endl;
422 #ifdef G4ANALYSIS_USE
423  if(fNtup[i]) fNtup[i]->addRow();
424 #endif
425 
426 #ifdef G4ANALYSIS_USE_ROOT
427  float *ar=new float[fRcol[i]];
428  for (G4int j=0; j < fRcol[i]; j++) {
429 // G4cout << i << " " << fRarray[i][j] << G4endl;
430  ar[j] = fRarray[i][j];
431  }
432  if(fROOTntup[i]) fROOTntup[i]->Fill(ar);
433  delete ar;
434 #endif
435 
436 }
437 
438 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
439 
441 {
442  fHistName = nam;
443 }
444 
445 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
446 
448 {
449  return fHistName;
450 }
451 
452 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
453 
455 {
456  fHistType = nam;
457 }
458 
459 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
460 
462 {
463  return fHistType;
464 }
465 
466 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
467 
std::vector< G4double > fXmin
Definition: exrdmHisto.hh:143
const G4String & GetFileName() const
Definition: exrdmHisto.cc:447
void FillTuple(G4int, const G4String &, G4double)
Definition: exrdmHisto.cc:375
G4String & remove(str_size)
Definition of the exrdmHisto class.
const G4String & FileType() const
Definition: exrdmHisto.cc:461
G4String name
Definition: TRTMaterials.hh:40
void SetFileType(const G4String &)
Definition: exrdmHisto.cc:454
std::vector< G4String > fTupleId
Definition: exrdmHisto.hh:149
void FillHisto(G4int, G4double, G4double)
Definition: exrdmHisto.cc:274
Definition of the exrdmHistoMessenger class.
void AddTuple(const G4String &, const G4String &, const G4String &)
Definition: exrdmHisto.cc:333
G4String fHistName
Definition: exrdmHisto.hh:117
std::vector< G4String > fTitles
Definition: exrdmHisto.hh:147
void SetHisto1D(G4int, G4int, G4double, G4double, G4double)
Definition: exrdmHisto.cc:254
G4String fHistType
Definition: exrdmHisto.hh:118
int G4int
Definition: G4Types.hh:78
void SetFileName(const G4String &)
Definition: exrdmHisto.cc:440
void AddRow(G4int)
Definition: exrdmHisto.cc:419
G4GLOB_DLL std::ostream G4cout
std::vector< G4int > fActive
Definition: exrdmHisto.hh:141
std::vector< G4String > fIds
Definition: exrdmHisto.hh:146
G4int fNTuple
Definition: exrdmHisto.hh:121
G4int fDefaultAct
Definition: exrdmHisto.hh:123
G4bool contains(const std::string &) const
std::vector< G4double > fUnit
Definition: exrdmHisto.hh:145
G4int fNHisto
Definition: exrdmHisto.hh:120
void Book()
Definition: exrdmHisto.cc:113
std::vector< G4double > fXmax
Definition: exrdmHisto.hh:144
std::vector< G4String > fTupleListROOT
Definition: exrdmHisto.hh:151
std::vector< G4String > fTupleName
Definition: exrdmHisto.hh:148
#define G4endl
Definition: G4ios.hh:61
std::vector< G4String > fTupleList
Definition: exrdmHisto.hh:150
double G4double
Definition: G4Types.hh:76
void ScaleHisto(G4int, G4double)
Definition: exrdmHisto.cc:301
G4int fVerbose
Definition: exrdmHisto.hh:122
void Add1D(const G4String &, const G4String &, G4int nb=100, G4double x1=0., G4double x2=1., G4double u=1.)
Definition: exrdmHisto.cc:226
void Save()
Definition: exrdmHisto.cc:196
exrdmHistoMessenger * fMessenger
Definition: exrdmHisto.hh:139
std::vector< G4int > fBins
Definition: exrdmHisto.hh:142