Geant4  10.03
plotG.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 #define USE_CANVASINTAB
27 
28 #ifdef USE_CANVASINTAB
29 #include "CanvasInTab.hh"
30 #endif
31 
32 #include <fstream>
33 #include <string>
34 #include <sstream>
35 #include <vector>
36 #include <set>
37 #include <map>
38 #include <iostream>
39 #include <iomanip>
40 #include <locale>
41 #include <sstream>
42 #include <cstring>
43 
44 #include <TApplication.h>
45 #include <TGApplication.h>
46 #include <TTree.h>
47 #include <TNtuple.h>
48 #include <TBranch.h>
49 #include <TProfile.h>
50 #include <TGraph.h>
51 #include <TAxis.h>
52 #include <TGraphErrors.h>
53 #include <TROOT.h>
54 #include <TCanvas.h>
55 #include <TFile.h>
56 #include <TGFileBrowser.h>
57 #include <TGFileDialog.h>
58 #include <TChain.h>
59 #include <TColor.h>
60 using namespace std;
61 
62 //------------------------------------------------------------------------------
63 
64 const TGFileInfo* OpenRootFile()
65 {
66  const char *gOpenAsTypes[] = {
67  "ROOT files", "*.root",
68  "All files", "*"
69  };
70 
71  static TGFileInfo fi;
72  fi.fFileTypes = gOpenAsTypes;
73  // fi.SetMultipleSelection(kTRUE);
74  // User must check the box "multiple selection" in the dialog box
75  // fi.fIniDir = StrDup(".");
76  new TGFileDialog(gClient->GetRoot(),
77  gClient->GetRoot(), kFDOpen, &fi);
78 
79  return &fi;
80 }
81 
82 //------------------------------------------------------------------------------
83 
85 {
87  {
88  fNEvent = 0;
89  fNumber = 0;
90  fG = 0.;
91  fG2 = 0.;
92  }
93 
94  SpeciesInfoAOS(const SpeciesInfoAOS& right) // Species A(B);
95  {
96  fNEvent = right.fNEvent;
97  fNumber = right.fNumber;
98  fG = right.fG;
99  fG2 = right.fG2;
100  fName = right.fName;
101  }
102 
104  {
105  if(&right == this) return *this;
106  fNEvent = right.fNEvent;
107  fNumber = right.fNumber;
108  fG = right.fG;
109  fG2 = right.fG2;
110  fName = right.fName;
111  return *this;
112  }
113 
114  int fNEvent;
115  int fNumber;
116  double fG;
117  double fG2;
118  string fName;
119 };
120 
121 //------------------------------------------------------------------------------
122 
124 {
126  {
127  fRelatErr = 0;
128  }
129 
131  fG(right.fG),
132  fGerr(right.fGerr),
133  fTime(right.fTime),
134  fRelatErr(right.fRelatErr),
135  fName(right.fName)
136  {}
137 
139  {
140  if(this == &right) return *this;
141  fG = right.fG;
142  fGerr = right.fGerr;
143  fTime = right.fTime;
144  fRelatErr = right.fRelatErr;
145  fName = right.fName;
146  return *this;
147  }
148 
149  std::vector<double> fG;
150  std::vector<double> fGerr;
151  std::vector<double> fTime;
152  double fRelatErr;
153  string fName;
154 };
155 
156 //------------------------------------------------------------------------------
157 
158 void ProcessSingleFile(TFile* file)
159 {
160  int speciesID;
161  int number;
162  int nEvent;
163  char speciesName[500];
164  double time; // time
165  double sumG; // sum of G over all events
166  double sumG2; // sum of G^2 over all events
167 
168  TTree* tree = (TTree*)file->Get("species");
169  tree->SetBranchAddress("speciesID", &speciesID);
170  tree->SetBranchAddress("number", &number);
171  tree->SetBranchAddress("nEvent", &nEvent);
172  tree->SetBranchAddress("speciesName", &speciesName);
173  tree->SetBranchAddress("time", &time);
174  tree->SetBranchAddress("sumG", &sumG);
175  tree->SetBranchAddress("sumG2", &sumG2);
176 
177  Long64_t nentries = tree->GetEntries();
178  // cout << nentries <<" entries" << endl;
179 
180  if(nentries == 0)
181  {
182  cout << "No entries found in the tree species contained in the file "
183  << file->GetPath() << endl;
184  exit(1);
185  }
186 
187  //----------------------------------------------------------------------------
188  // This first loop is used in case the processed ROOT file is issued from the
189  // accumulation of several ROOT files (e.g. hadd)
190 
191  std::map<int, std::map<double, SpeciesInfoAOS>> speciesTimeInfo;
192 
193  for (int j=0; j < nentries; j++)
194  {
195  tree->GetEntry(j);
196 
197  SpeciesInfoAOS& infoAOS = speciesTimeInfo[speciesID][time];
198 
199  infoAOS.fNumber += number;
200  infoAOS.fG += sumG;
201  infoAOS.fG2 += sumG2;
202  infoAOS.fNEvent += nEvent;
203  infoAOS.fName = speciesName;
204  }
205 
206  //----------------------------------------------------------------------------
207 
208  std::map<int, SpeciesInfoSOA> speciesInfo;
209 
210  auto it_SOA = speciesTimeInfo.begin();
211  auto end_SOA = speciesTimeInfo.end();
212 
213  for (; it_SOA!=end_SOA ; ++it_SOA)
214  {
215  const int _speciesID = it_SOA->first;
216  SpeciesInfoSOA& info = speciesInfo[_speciesID];
217 
218  auto it2 = it_SOA->second.begin();
219  auto end2 = it_SOA->second.end();
220 
221  info.fName = it2->second.fName;
222  const size_t size2 = it_SOA->second.size();
223  info.fG.resize(size2);
224  info.fGerr.resize(size2);
225  info.fTime.resize(size2);
226 
227  for(int i2 = 0 ;it2!=end2;++it2, ++i2)
228  {
229  SpeciesInfoAOS& infoAOS = it2->second;
230 
231  double _SumG2 = infoAOS.fG2;
232  double _MeanG = infoAOS.fG/infoAOS.fNEvent;
233  double _Gerr = sqrt((_SumG2/infoAOS.fNEvent - pow(_MeanG,2))
234  /(infoAOS.fNEvent-1) );
235 
236  info.fG[i2] = _MeanG;
237  info.fGerr[i2] = _Gerr;
238  info.fTime[i2] = it2->first;
239  info.fRelatErr += _Gerr/(_MeanG + 1e-30); // add an epsilon to prevent NAN
240  }
241  }
242 
243  //----------------------------------------------------------------------------
244 
245 #ifdef USE_CANVASINTAB
246  CanvasInTab* myFrame =
247  new CanvasInTab(gClient->GetRoot(), 500, 500);
248 #endif
249 
250  std::map<int, SpeciesInfoSOA>::iterator it = speciesInfo.begin();
251  std::map<int, SpeciesInfoSOA>::iterator end = speciesInfo.end();
252 
253  for (; it != end; ++it)
254  {
255  speciesID = it->first;
256  SpeciesInfoSOA& info = it->second;
257 // if(strstr(info.fName.c_str(), "H2O^") != 0) continue;
258 
259  if(info.fG.empty()) continue;
260 
261  TGraphErrors* gSpecies = new TGraphErrors(info.fG.size(),
262  info.fTime.data(),
263  info.fG.data(),
264  0,
265  info.fGerr.data());
266 
267 #ifdef USE_CANVASINTAB
268  int nCanvas = myFrame->AddCanvas(info.fName.c_str());
269  myFrame->GetCanvas(nCanvas);
270  TCanvas* cSpecies = myFrame->GetCanvas(nCanvas);
271 #else
272  TCanvas* cSpecies = new TCanvas(info.fName.c_str(),
273  info.fName.c_str());
274 #endif
275 
276  cSpecies->cd();
277  int color = (2+speciesID)%TColor::GetNumberOfColors();
278  if(color == 5 || color==10 || color==0) ++color;
279 
280  // cout << info.fName.c_str() << " " << color << endl;
281 
282  gSpecies->SetMarkerStyle(20+speciesID);
283  gSpecies->SetMarkerColor(color);
284  info.fRelatErr /= (double)info.fG.size();
285 
286  gSpecies->SetTitle((info.fName
287  + " - speciesID: "
288  + std::to_string(speciesID)+" rel. Err. "
289  + std::to_string(info.fRelatErr)).c_str() );
290  gSpecies->GetXaxis()->SetTitle("Time [ns]");
291  gSpecies->GetYaxis()->SetTitle("G [molecules/100 eV]");
292  gSpecies->Draw("ap");
293  cSpecies->SetLogx();
294  }
295 
296 #ifdef USE_CANVASINTAB
297  int nCanvas = myFrame->GetNCanvas();
298  for(int i = 0 ; i < nCanvas ; ++i)
299  {
300  myFrame->GetCanvas(i)->Update();
301  }
302 #endif
303 }
304 
305 //------------------------------------------------------------------------------
306 
307 int ProcessSingleFile(const char* filePath)
308 {
309  if(filePath == 0 || strlen(filePath) == 0)
310  {
311  perror("You must provide a valid file");
312  return 1;
313  }
314 
315  TFile* file = TFile::Open(filePath);
316 
317  if(file == 0)
318  {
319  perror ("Error opening ntuple file");
320  exit(1);
321  }
322 
323  if(!file-> IsOpen())
324  {
325  perror ("Error opening ntuple file");
326  exit(1);
327  }
328  else
329  {
330  cout << "Opening ntple file " << filePath << endl;
331  }
332  ProcessSingleFile(file);
333  return 0;
334 }
335 
336 //------------------------------------------------------------------------------
337 
338 #define _PROCESS_ONE_FILE_ ProcessSingleFile
339 //#define _PROCESS_ONE_FILE_ ProcessSingleFileTProfile
340 
341 int main(int argc, char **argv)
342 {
343  //--------------------------------
344  int initialArgc = argc;
345  vector<char*> initialArgv(argc);
346  for(int i = 0 ; i < argc ; ++i)
347  {
348  initialArgv[i] = argv[i];
349  }
350  //--------------------------------
351 
352  TApplication* rootApp = new TApplication("PlotG",&argc, argv);
353 
354  const char* filePath = 0;
355 
356  if(initialArgc == 1) // no file provided in argument
357  {
358  const TGFileInfo* fileInfo = OpenRootFile();
359  filePath = fileInfo->fFilename;
360  if(fileInfo->fFileNamesList && fileInfo->fFileNamesList->GetSize()>1)
361  {
362  // several files selected
363  // user has to tick "Multiple selection"
364  perror("Multiple selection of files not supported, implement your own!");
365  //
366  // For instance, start from:
367  // TChain* tree = new TChain("species");
368  // tree->AddFileInfoList(fileInfo->fFileNamesList);
369  // Or call ProcessSingleFile for each file,
370  // you'll need to do some adaptation
371  }
372  else
373  {
374  if(_PROCESS_ONE_FILE_(filePath)) return 1;
375  }
376  }
377  else // a file is provided in argument
378  {
379  filePath = initialArgv[1];
380  if(_PROCESS_ONE_FILE_(filePath)) return 1;
381  }
382 
383  rootApp->Run();
384  delete rootApp;
385  return 0;
386 }
void ProcessSingleFile(TFile *file)
Definition: plotG.cc:158
SpeciesInfoSOA(const SpeciesInfoSOA &right)
Definition: plotG.cc:130
std::vector< double > fTime
Definition: plotG.cc:151
SpeciesInfoAOS()
Definition: plotG.cc:86
TCanvas * GetCanvas(int i)
Definition: CanvasInTab.cc:140
string fName
Definition: plotG.cc:153
G4String fName
Definition: G4AttUtils.hh:55
double fRelatErr
Definition: plotG.cc:152
int main(int argc, char **argv)
Definition: plotG.cc:341
SpeciesInfoSOA()
Definition: plotG.cc:125
string fName
Definition: plotG.cc:118
std::vector< double > fGerr
Definition: plotG.cc:150
size_t GetNCanvas() const
Definition: CanvasInTab.hh:49
#define _PROCESS_ONE_FILE_
Definition: plotG.cc:338
size_t AddCanvas(const char *name="New tab")
Definition: CanvasInTab.cc:121
SpeciesInfoAOS & operator=(const SpeciesInfoAOS &right)
Definition: plotG.cc:103
const TGFileInfo * OpenRootFile()
Definition: plotG.cc:64
SpeciesInfoAOS(const SpeciesInfoAOS &right)
Definition: plotG.cc:94
double fG2
Definition: plotG.cc:117
std::vector< double > fG
Definition: plotG.cc:149
SpeciesInfoSOA & operator=(const SpeciesInfoSOA &right)
Definition: plotG.cc:138
double fG
Definition: plotG.cc:116