| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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",
65
66 static TGFileInfo fi;
67 fi.fFileTypes = gOpenAsTypes;
68 // fi.SetMultipleSelection(kTRUE);
69 // User must check the box "multiple selecti
70 // fi.fIniDir = StrDup(".");
71 new TGFileDialog(gClient->GetRoot(), gClient
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
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 SpeciesInf
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 SpeciesInf
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", &species
161 tree->SetBranchAddress("number", &number);
162 tree->SetBranchAddress("nEvent", &nEvent);
163 tree->SetBranchAddress("speciesName", &speci
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 spec
173 << endl;
174 exit(1);
175 }
176
177 //------------------------------------------
178 // This first loop is used in case the proce
179 // accumulation of several ROOT files (e.g.
180
181 std::map<int, std::map<double, SpeciesInfoAO
182
183 for (int j = 0; j < nentries; j++) {
184 tree->GetEntry(j);
185
186 SpeciesInfoAOS& infoAOS = speciesTimeInfo[
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[_specie
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.fNE
220 double _Gerr = sqrt((_SumG2 / infoAOS.fN
221
222 info.fG[i2] = _MeanG;
223 info.fGerr[i2] = _Gerr;
224 info.fTime[i2] = it2->first;
225 info.fRelatErr += _Gerr / (_MeanG + 1e-3
226 }
227 }
228
229 //------------------------------------------
230
231 #ifdef USE_CANVASINTAB
232 CanvasInTab* myFrame = new CanvasInTab(gClie
233 #endif
234
235 std::map<int, SpeciesInfoSOA>::iterator it =
236 std::map<int, SpeciesInfoSOA>::iterator end
237
238 for (; it != end; ++it) {
239 speciesID = it->first;
240 SpeciesInfoSOA& info = it->second;
241 // if(strstr(info.fName.c_str(), "H2O^"
242
243 if (info.fG.empty()) continue;
244
245 TGraphErrors* gSpecies =
246 new TGraphErrors(info.fG.size(), info.fT
247
248 #ifdef USE_CANVASINTAB
249 int nCanvas = myFrame->AddCanvas(info.fNam
250 myFrame->GetCanvas(nCanvas);
251 TCanvas* cSpecies = myFrame->GetCanvas(nCa
252 #else
253 TCanvas* cSpecies = new TCanvas(info.fName
254 #endif
255
256 cSpecies->cd();
257 int color = (2 + speciesID) % TColor::GetN
258 if (color == 5 || color == 10 || color ==
259
260 // cout << info.fName.c_str() << " " << co
261
262 gSpecies->SetMarkerStyle(20 + speciesID);
263 gSpecies->SetMarkerColor(color);
264 info.fRelatErr /= (double)info.fG.size();
265
266 gSpecies->SetTitle((info.fName + " - speci
267 + std::to_string(info.
268 .c_str());
269 gSpecies->GetXaxis()->SetTitle("Time [ns]"
270 gSpecies->GetYaxis()->SetTitle("G [molecul
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
305 }
306 ProcessSingleFile(file);
307 return 0;
308 }
309
310 //--------------------------------------------
311
312 #define _PROCESS_ONE_FILE_ ProcessSingleFile
313 // #define _PROCESS_ONE_FILE_ ProcessSingleFil
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("Pl
326
327 const char* filePath = 0;
328
329 if (initialArgc == 1) // no file provided i
330 {
331 const TGFileInfo* fileInfo = OpenRootFile(
332 filePath = fileInfo->fFilename;
333 if (fileInfo->fFileNamesList && fileInfo->
334 // several files selected
335 // user has to tick "Multiple selection"
336 perror("Multiple selection of files not
337 //
338 // For instance, start from:
339 // TChain* tree = new TChain("species"
340 // tree->AddFileInfoList(fileInfo->fFi
341 // Or call ProcessSingleFile for each fi
342 // you'll need to do some adaptation
343 }
344 else {
345 if (_PROCESS_ONE_FILE_(filePath)) return
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