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