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