src/erhic/Forester.cxx

0001 /**
0002  \file
0003  Implementation of class erhic::Forester.
0004
0005  \author    Thomas Burton
0006  \date      2011-06-23
0007  \copyright 2011 Brookhaven National Lab
0008  */
0009
0010 #include "eicsmear/erhic/Forester.h"
0011
0012 #include <iomanip>
0013 #include <memory>
0014 #include <stdexcept>
0015 #include <string>
0016
0017 #include <TRefArray.h>
0018 #include <TString.h>
0019
0020 #include "eicsmear/erhic/EventFactory.h"
0021 #include "eicsmear/erhic/File.h"
0022 #include "eicsmear/erhic/ParticleIdentifier.h"
0023
0024 #include "eicsmear/gzstream.h"
0025
0026 namespace erhic {
0027
0028 Forester::Forester()
0029 : mQuit(false)
0030 , mVerbose(false)
0031 , mTree(NULL)
0032 , mEvent(NULL)
0033 , mFile(NULL)
0034 , mRootFile(NULL)
0035 , mMaxNEvents(0)
0036 , mInterval(1)
0037 , mTextFile(NULL)
0038 , mInputName("default.txt")
0039 , mOutputName("default.root")
0040 , mTreeName("EICTree")
0041 , mBranchName("event")
0042 , mFactory(NULL) {
0043 }
0044
0045 Forester::~Forester() {
0046   if (mFile) {
0047     delete mFile;
0048     mFile = NULL;
0049   }  // if
0050   if (mEvent) {
0051     delete mEvent;
0052     mEvent = NULL;
0053   }  // if
0054   if (mFactory) {
0055     delete mFactory;
0056     mFactory = NULL;
0057   }  // if
0058   if (mRootFile) {
0059     delete mRootFile;
0060     mRootFile = NULL;
0061   }  // if
0062   // if (mTextFile) {
0063   //   delete mTextFile;
0064   //   mTextFile = NULL;
0065   // }  // if
0066   // We don't delete the mTree pointer because mRootFile
0067   // has ownership of it.
0068 }
0069
0070 Long64_t Forester::Plant() {
0071   try {
0072     // Initialisation of the input and output files.
0073     OpenInput();
0074     SetupOutput();
0075     if (BeVerbose()) {
0076       std::cout << "\nProcessing " << GetInputFileName() << std::endl;
0077     }  // if
0078     /**
0079      \todo Get rid of the static counter. Replace with a member
0080      that is reset every time Plant() is called.
0081      */
0082     static int i(0);
0083     while (!MustQuit()) {
0084       ++i;
0085       if (BeVerbose() && i % mInterval == 0) {
0086         // Make the field just wide enough for the maximum
0087         // number of events.
0088         int width = static_cast<int>(::log10(GetMaxNEvents()) + 1);
0089         std::cout << "Processing event "<< std::setw(width) << i;
0090         if (GetMaxNEvents() > 0) {
0091           std::cout << "/" << std::setw(width) << GetMaxNEvents();
0092         }  // if
0093         std::cout << std::endl;
0094       }  // if
0095       // Build the next event
0096       if (mEvent) {
0097         delete mEvent;
0098         mEvent = NULL;
0099       }  // if
0100       // Catch exceptions from event builder here so we don't break
0101       // out of the whole tree building loop for a single bad event.
0102       try {
0103         mEvent = mFactory->Create();
0104         // Fill the tree
0105         if (mEvent) {
0106           mTree->Fill();
0107           if (GetMaxNEvents() > 0 && i >= GetMaxNEvents()) {
0108             SetMustQuit(true);  // Hit max number of events, so quit
0109           }  // if
0110           mStatus.ModifyEventCount(1);
0111           mStatus.ModifyParticleCount(mEvent->GetNTracks());
0112           // We must ResetBranchAddress before deleting the event.
0113         } else {
0114           break;
0115         }  // if
0116       }  // try
0117       catch(std::exception& e) {
0118         std::cerr << "Caught exception in Forester::Plant(): "
0119         << e.what() << std::endl;
0120         std::cerr << "Event will be skipped..." << std::endl;
0121       }  // catch
0122     }  // while
0123     Finish();
0124     return 0;
0125   }  // try
0126   catch(std::exception& e) {
0127     std::cerr << "Caught exception in Forester::Plant(): "
0128     << e.what() << std::endl;
0129     return -1;
0130   }  // catch
0131 }
0132
0133 bool Forester::OpenInput() {
0134   try {
0135     // Open the input file for reading.
0136     if ( TString(GetInputFileName()).EndsWith("gz", TString::kIgnoreCase) ||
0137      TString(GetInputFileName()).EndsWith("zip", TString::kIgnoreCase)){
0138       auto tmp = std::make_shared<igzstream>();
0139       tmp->open(GetInputFileName().c_str());
0140       // Throw a runtime_error if the file could not be opened.
0141       if (!tmp->good()) {
0142     std::string message("Unable to open file ");
0143     throw std::runtime_error(message.append(GetInputFileName()));
0144       }  // if
0145       mTextFile = tmp;
0146     } else {
0147       auto tmp = std::make_shared<std::ifstream>();
0148       tmp->open(GetInputFileName().c_str());
0149       // Throw a runtime_error if the file could not be opened.
0150       if (!tmp->good()) {
0151     std::string message("Unable to open file ");
0152     throw std::runtime_error(message.append(GetInputFileName()));
0153       }  // if
0154       mTextFile = tmp;
0155     }
0156
0157     // Determine which Monte Carlo generator produced the file.
0158     // hand over file name too, because the next function is destructive to the stream
0159     // and gzipped hepmc files need to reopen the stream
0160     mFile = erhic::FileFactory::GetInstance().GetFile(mTextFile, GetInputFileName());
0161     if (!mFile) {
0162       throw std::runtime_error(GetInputFileName() +
0163                        " is not from a supported generator");
0164     }  // for
0165     mFactory = mFile->CreateEventFactory(*mTextFile);
0166     mFactory->mAdditionalInformation["generator"]=mFile->GetGeneratorName();
0167     return true;
0168   } // Pass the exception on to be dealt with higher up the food chain.
0169   catch(std::exception&) {
0170     throw;
0171   }  // catch
0172 }
0173
0174 bool Forester::SetupOutput() {
0175   try {
0176     // Open the ROOT file and check it opened OK
0177     mRootFile = new TFile(GetOutputFileName().c_str(), "RECREATE");
0178     if (!mRootFile->IsOpen()) {
0179       std::string message("Unable to open file ");
0180       throw std::runtime_error(message.append(GetOutputFileName()));
0181     }  // if
0182     // Create the tree and check for errors
0183     mTree = new TTree(GetTreeName().c_str(), "my EIC tree");
0184     if (!mTree) {
0185       std::string message("Error allocating TTree ");
0186       throw std::runtime_error(message.append(GetTreeName()));
0187     }  // if
0188     // Allocate memory for the branch buffer and
0189     // add the branch to the tree
0190     AllocateEvent();
0191     mTree->Branch(GetBranchName().c_str(), mEvent->ClassName(),
0192                   &mEvent, 32000, 99);
0193     // Auto-save every 500 MB
0194     mTree->SetAutoSave(500LL * 1024LL * 1024LL);
0195     // Align the input file at the start of the first event (event generator dependent).
0196     mFactory->FindFirstEvent();
0197     // Start timing after opening and creating files,
0198     // before looping over events
0199     mStatus.StartTimer();
0200     return true;
0201   }  // try...
0202   catch(std::exception&) {
0203     throw;
0204   }  // catch
0205 }
0206
0207 void Forester::Finish() {
0208   if (BeVerbose()) {
0209     std::cout << "\nProcessed " << GetInputFileName() << std::endl;
0210   }  // if
0211   // Write the TTree to the file.
0212   mRootFile = mTree->GetCurrentFile();
0213   mRootFile->cd();
0214   mTree->Write();
0215   for ( auto namedobject : mFactory->mObjectsToWriteAtTheEnd ){
0216     namedobject.second->Write(namedobject.first);
0217   }
0218   mRootFile->ls();
0219   // Write the Forester itself to make it easier to reproduce the file
0220   // with the same settings.
0221   Write("forester");
0222   // Reset quit flag in case of further runs.
0223   SetMustQuit(false);
0224   // Stop timing the run.
0225   mStatus.StopTimer();
0226   if (BeVerbose()) {
0227     GetGetStatus().Print(std::cout);  // Messages for the user
0228   }  // if
0229   mRootFile->Close();
0230 }
0231
0232 bool Forester::AllocateEvent() {
0233   try {
0234     if (mEvent) {
0235       delete mEvent;
0236       mEvent = NULL;
0237     }  // if
0238     mEvent = mFile->AllocateEvent();
0239     return mEvent;
0240   }  // try...
0241   // Catch exceptions and pass up the food chain
0242   catch(std::exception&) {
0243     throw;
0244   }  // catch
0245 }
0246
0247 bool Forester::FindFirstEvent() {
0248   // Naughty kludge alert!
0249   // The first line was already read to determine the generator.
0250   // The header in the text files is six lines, so
0251   // read the remaining five lines of header.
0252   std::getline(*mTextFile, mLine);
0253   std::getline(*mTextFile, mLine);
0254   std::getline(*mTextFile, mLine);
0255   std::getline(*mTextFile, mLine);
0256   std::getline(*mTextFile, mLine);
0257   return true;
0258 }
0259
0260 void Forester::Print(std::ostream& os) const {
0261   os << "Input file: " << mInputName << std::endl;
0262   os << "Output file: " << mOutputName << std::endl;
0263   os << "Output tree: " << mTreeName << std::endl;
0264   os << "Output branch: " << mBranchName << std::endl;
0265   os << "Maximum number of events: " << mMaxNEvents << std::endl;
0266   if (mEvent) {
0267     os << "Event type: " << mEvent->ClassName() << std::endl;
0268   }  // if
0269 }
0270
0271 void Forester::Print(Option_t* /* not used */) const {
0272   Print(std::cout);
0273 }
0274
0275   Forester::Status::Status()
0276     : mNEvents(0)
0277     , mNParticles(0) {
0278     // Initialise the start and end time to the creation time and reset
0279     // the timer to ensure it is at zero.
0280     std::time(&mStartTime);
0281     mEndTime = mStartTime;
0282     mTimer.Reset();
0283   }
0284
0285   Forester::Status::~Status() { /* noop */ }
0286
0287   std::ostream& Forester::Status::Print(std::ostream& os) const {
0288     // Put start and end times in different os <<... otherwise I get
0289     // the same time for each...
0290     os << "Began on " << std::ctime(&mStartTime);
0291     os << "Ended on " << std::ctime(&mEndTime);
0292     os << "Processed " << mNEvents << " events containing "
0293        << mNParticles << " particles in "
0294        << mTimer.RealTime() << " seconds "
0295        << '(' << mTimer.RealTime()/mNEvents <<" sec/event)" << std::endl;
0296     return os;
0297   }
0298
0299   void Forester::Status::StartTimer() {
0300     std::time(&mStartTime);
0301     mTimer.Start();
0302   }
0303
0304   void Forester::Status::StopTimer() {
0305     std::time(&mEndTime);
0306     mTimer.Stop();
0307   }
0308
0309   void Forester::Status::ModifyEventCount(Long64_t count) {
0310     mNEvents += count;
0311   }
0312
0313   void Forester::Status::ModifyParticleCount(Long64_t count) {
0314     mNParticles += count;
0315   }
0316
0317   // ClassImp( ForesterStatus ); // throws error for some reason
0318
0319
0320
0321 }  // namespace erhic