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