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