chem4/plot/plotG.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /// \file plotG.cc
0027 /// \brief Main program of the chem4/plot example
0028
0029 #define USE_CANVASINTAB
0030
0031 #ifdef USE_CANVASINTAB
0032 #  include "CanvasInTab.hh"
0033 #endif
0034
0035 #include <TApplication.h>
0036 #include <TAxis.h>
0037 #include <TBranch.h>
0038 #include <TCanvas.h>
0039 #include <TChain.h>
0040 #include <TColor.h>
0041 #include <TFile.h>
0042 #include <TGApplication.h>
0043 #include <TGFileBrowser.h>
0044 #include <TGFileDialog.h>
0045 #include <TGraph.h>
0046 #include <TGraphErrors.h>
0047 #include <TNtuple.h>
0048 #include <TProfile.h>
0049 #include <TROOT.h>
0050 #include <TTree.h>
0051 #include <cstring>
0052 #include <fstream>
0053 #include <iomanip>
0054 #include <iostream>
0055 #include <locale>
0056 #include <map>
0057 #include <set>
0058 #include <sstream>
0059 #include <string>
0060 #include <vector>
0061 using namespace std;
0062
0063 //------------------------------------------------------------------------------
0064
0065 const TGFileInfo* OpenRootFile()
0066 {
0067   const char* gOpenAsTypes[] = {"ROOT files", "*.root", "All files", "*"};
0068
0069   static TGFileInfo fi;
0070   fi.fFileTypes = gOpenAsTypes;
0071   //  fi.SetMultipleSelection(kTRUE);
0072   // User must check the box "multiple selection" in the dialog box
0073   //  fi.fIniDir = StrDup(".");
0074   new TGFileDialog(gClient->GetRoot(), gClient->GetRoot(), kFDOpen, &fi);
0075
0076   return &fi;
0077 }
0078
0079 //------------------------------------------------------------------------------
0080
0081 struct SpeciesInfoAOS
0082 {
0083     SpeciesInfoAOS()
0084     {
0085       fNEvent = 0;
0086       fNumber = 0;
0087       fG = 0.;
0088       fG2 = 0.;
0089     }
0090
0091     SpeciesInfoAOS(const SpeciesInfoAOS& right)  // Species A(B);
0092     {
0093       fNEvent = right.fNEvent;
0094       fNumber = right.fNumber;
0095       fG = right.fG;
0096       fG2 = right.fG2;
0097       fName = right.fName;
0098     }
0099
0100     SpeciesInfoAOS& operator=(const SpeciesInfoAOS& right)  // A = B
0101     {
0102       if (&right == this) return *this;
0103       fNEvent = right.fNEvent;
0104       fNumber = right.fNumber;
0105       fG = right.fG;
0106       fG2 = right.fG2;
0107       fName = right.fName;
0108       return *this;
0109     }
0110
0111     int fNEvent;
0112     int fNumber;
0113     double fG;
0114     double fG2;
0115     string fName;
0116 };
0117
0118 //------------------------------------------------------------------------------
0119
0120 struct SpeciesInfoSOA
0121 {
0122     SpeciesInfoSOA() { fRelatErr = 0; }
0123
0124     SpeciesInfoSOA(const SpeciesInfoSOA& right)
0125       : fG(right.fG),
0126         fGerr(right.fGerr),
0127         fTime(right.fTime),
0128         fRelatErr(right.fRelatErr),
0129         fName(right.fName)
0130     {}
0131
0132     SpeciesInfoSOA& operator=(const SpeciesInfoSOA& right)
0133     {
0134       if (this == &right) return *this;
0135       fG = right.fG;
0136       fGerr = right.fGerr;
0137       fTime = right.fTime;
0138       fRelatErr = right.fRelatErr;
0139       fName = right.fName;
0140       return *this;
0141     }
0142
0143     std::vector<double> fG;
0144     std::vector<double> fGerr;
0145     std::vector<double> fTime;
0146     double fRelatErr;
0147     string fName;
0148 };
0149
0150 //------------------------------------------------------------------------------
0151
0152 void ProcessSingleFile(TFile* file)
0153 {
0154   int speciesID;
0155   int number;
0156   int nEvent;
0157   char speciesName[500];
0158   double time;  // time
0159   double sumG;  // sum of G over all events
0160   double sumG2;  // sum of G^2 over all events
0161
0162   TTree* tree = (TTree*)file->Get("species");
0163   tree->SetBranchAddress("speciesID", &speciesID);
0164   tree->SetBranchAddress("number", &number);
0165   tree->SetBranchAddress("nEvent", &nEvent);
0166   tree->SetBranchAddress("speciesName", &speciesName);
0167   tree->SetBranchAddress("time", &time);
0168   tree->SetBranchAddress("sumG", &sumG);
0169   tree->SetBranchAddress("sumG2", &sumG2);
0170
0171   Long64_t nentries = tree->GetEntries();
0172   // cout << nentries <<" entries" << endl;
0173
0174   if (nentries == 0) {
0175     cout << "No entries found in the tree species contained in the file " << file->GetPath()
0176          << endl;
0177     exit(1);
0178   }
0179
0180   //----------------------------------------------------------------------------
0181   // This first loop is used in case the processed ROOT file is issued from the
0182   // accumulation of several ROOT files (e.g. hadd)
0183
0184   std::map<int, std::map<double, SpeciesInfoAOS>> speciesTimeInfo;
0185
0186   for (int j = 0; j < nentries; j++) {
0187     tree->GetEntry(j);
0188
0189     SpeciesInfoAOS& infoAOS = speciesTimeInfo[speciesID][time];
0190
0191     infoAOS.fNumber += number;
0192     infoAOS.fG += sumG;
0193     infoAOS.fG2 += sumG2;
0194     infoAOS.fNEvent += nEvent;
0195     infoAOS.fName = speciesName;
0196   }
0197
0198   //----------------------------------------------------------------------------
0199
0200   std::map<int, SpeciesInfoSOA> speciesInfo;
0201
0202   auto it_SOA = speciesTimeInfo.begin();
0203   auto end_SOA = speciesTimeInfo.end();
0204
0205   for (; it_SOA != end_SOA; ++it_SOA) {
0206     const int _speciesID = it_SOA->first;
0207     SpeciesInfoSOA& info = speciesInfo[_speciesID];
0208
0209     auto it2 = it_SOA->second.begin();
0210     auto end2 = it_SOA->second.end();
0211
0212     info.fName = it2->second.fName;
0213     const size_t size2 = it_SOA->second.size();
0214     info.fG.resize(size2);
0215     info.fGerr.resize(size2);
0216     info.fTime.resize(size2);
0217
0218     for (int i2 = 0; it2 != end2; ++it2, ++i2) {
0219       SpeciesInfoAOS& infoAOS = it2->second;
0220
0221       double _SumG2 = infoAOS.fG2;
0222       double _MeanG = infoAOS.fG / infoAOS.fNEvent;
0223       double _Gerr = sqrt((_SumG2 / infoAOS.fNEvent - pow(_MeanG, 2)) / (infoAOS.fNEvent - 1));
0224
0225       info.fG[i2] = _MeanG;
0226       info.fGerr[i2] = _Gerr;
0227       info.fTime[i2] = it2->first;
0228       info.fRelatErr += _Gerr / (_MeanG + 1e-30);  // add an epsilon to prevent NAN
0229     }
0230   }
0231
0232   //----------------------------------------------------------------------------
0233
0234 #ifdef USE_CANVASINTAB
0235   CanvasInTab* myFrame = new CanvasInTab(gClient->GetRoot(), 500, 500);
0236 #endif
0237
0238   std::map<int, SpeciesInfoSOA>::iterator it = speciesInfo.begin();
0239   std::map<int, SpeciesInfoSOA>::iterator end = speciesInfo.end();
0240
0241   for (; it != end; ++it) {
0242     speciesID = it->first;
0243     SpeciesInfoSOA& info = it->second;
0244     //    if(strstr(info.fName.c_str(), "H2O^") != 0) continue;
0245
0246     if (info.fG.empty()) continue;
0247
0248     TGraphErrors* gSpecies =
0249       new TGraphErrors(info.fG.size(), info.fTime.data(), info.fG.data(), 0, info.fGerr.data());
0250
0251 #ifdef USE_CANVASINTAB
0252     int nCanvas = myFrame->AddCanvas(info.fName.c_str());
0253     myFrame->GetCanvas(nCanvas);
0254     TCanvas* cSpecies = myFrame->GetCanvas(nCanvas);
0255 #else
0256     TCanvas* cSpecies = new TCanvas(info.fName.c_str(), info.fName.c_str());
0257 #endif
0258
0259     cSpecies->cd();
0260     int color = (2 + speciesID) % TColor::GetNumberOfColors();
0261     if (color == 5 || color == 10 || color == 0) ++color;
0262
0263     // cout << info.fName.c_str() << " " << color << endl;
0264
0265     gSpecies->SetMarkerStyle(20 + speciesID);
0266     gSpecies->SetMarkerColor(color);
0267     info.fRelatErr /= (double)info.fG.size();
0268
0269     gSpecies->SetTitle((info.fName + " - speciesID: " + std::to_string(speciesID) + " rel. Err. "
0270                         + std::to_string(info.fRelatErr))
0271                          .c_str());
0272     gSpecies->GetXaxis()->SetTitle("Time [ns]");
0273     gSpecies->GetYaxis()->SetTitle("G [molecules/100 eV]");
0274     gSpecies->Draw("ap");
0275     cSpecies->SetLogx();
0276   }
0277
0278 #ifdef USE_CANVASINTAB
0279   int nCanvas = myFrame->GetNCanvas();
0280   for (int i = 0; i < nCanvas; ++i) {
0281     myFrame->GetCanvas(i)->Update();
0282   }
0283 #endif
0284 }
0285
0286 //------------------------------------------------------------------------------
0287
0288 int ProcessSingleFile(const char* filePath)
0289 {
0290   if (filePath == 0 || strlen(filePath) == 0) {
0291     perror("You must provide a valid file");
0292     return 1;
0293   }
0294
0295   TFile* file = TFile::Open(filePath);
0296
0297   if (file == 0) {
0298     perror("Error opening ntuple file");
0299     exit(1);
0300   }
0301
0302   if (!file->IsOpen()) {
0303     perror("Error opening ntuple file");
0304     exit(1);
0305   }
0306   else {
0307     cout << "Opening ntple file " << filePath << endl;
0308   }
0309   ProcessSingleFile(file);
0310   return 0;
0311 }
0312
0313 //------------------------------------------------------------------------------
0314
0315 #define _PROCESS_ONE_FILE_ ProcessSingleFile
0316 // #define _PROCESS_ONE_FILE_ ProcessSingleFileTProfile
0317
0318 int main(int argc, char** argv)
0319 {
0320   //--------------------------------
0321   int initialArgc = argc;
0322   vector<char*> initialArgv(argc);
0323   for (int i = 0; i < argc; ++i) {
0324     initialArgv[i] = argv[i];
0325   }
0326   //--------------------------------
0327
0328   TApplication* rootApp = new TApplication("PlotG", &argc, argv);
0329
0330   const char* filePath = 0;
0331
0332   if (initialArgc == 1)  // no file provided in argument
0333   {
0334     const TGFileInfo* fileInfo = OpenRootFile();
0335     filePath = fileInfo->fFilename;
0336     if (fileInfo->fFileNamesList && fileInfo->fFileNamesList->GetSize() > 1) {
0337       // several files selected
0338       // user has to tick "Multiple selection"
0339       perror("Multiple selection of files not supported, implement your own!");
0340       //
0341       // For instance, start from:
0342       //   TChain* tree = new TChain("species");
0343       //   tree->AddFileInfoList(fileInfo->fFileNamesList);
0344       // Or call ProcessSingleFile for each file,
0345       // you'll need to do some adaptation
0346     }
0347     else {
0348       if (_PROCESS_ONE_FILE_(filePath)) return 1;
0349     }
0350   }
0351   else  // a file is provided in argument
0352   {
0353     filePath = initialArgv[1];
0354     if (_PROCESS_ONE_FILE_(filePath)) return 1;
0355   }
0356
0357   rootApp->Run();
0358   delete rootApp;
0359   return 0;
0360 }