src/erhic/Pythia6.cxx

0001 /**
0002  \file
0003  Implementation of class erhic::Pythia6.
0004
0005  \author    Thomas Burton
0006  \date      2012-01-18
0007  \copyright 2012 Brookhaven National Lab
0008  */
0009
0010 #include "eicsmear/erhic/Pythia6.h"
0011
0012 #include <iostream>
0013 #include <stdexcept>
0014 #include <sstream>
0015 #include <string>
0016
0017 #include <TFile.h>
0018 #include <TObjString.h>
0019 #include <TPythia6.h>
0020 #include <TStopwatch.h>
0021 #include <TTree.h>
0022
0023 #include "eicsmear/erhic/EventMCFilterABC.h"
0024 #include "eicsmear/erhic/ParticleMC.h"
0025 #include "eicsmear/erhic/EventFactory.h"
0026
0027 namespace erhic {
0028
0029 Pythia6::Pythia6(TFile* file,
0030                  VirtualEventFactory* factory,
0031                  int nEvents,
0032                  const std::string& treeName,
0033                  const std::string& branchName,
0034                  int printInterval)
0035 : mPrintInterval(printInterval)
0036 , mFile(file)
0037 , mTree(NULL)
0038 , mNEvents(nEvents)
0039 , mNGenerated(0)
0040 , mNTrials(0)
0041 , mFactory(factory) {
0042   try {
0043     if (!file) {
0044       throw std::runtime_error("No file provided");
0045     }  // if
0046     if (!file->IsWritable()) {
0047       throw std::runtime_error("File is not writable");
0048     }  // if
0049     file->cd();
0050     mTree = new TTree(treeName.c_str(), "PYTHIA 6 events");
0051     mFactory->Branch(*mTree, branchName);
0052   }  // try
0053   catch(std::exception& e) {
0054     std::cout << "Caught exception in erhic::Pythia6::Pythia6() - " <<
0055     e.what() << std::endl;
0056     exit(1);
0057   }  // catch
0058 }
0059
0060 Pythia6::~Pythia6() {
0061   // Don't delete the file as that was externally provided
0062   // Don't delete the tree as that is associated with the file
0063   // and will be deleted when the file closes.
0064 }
0065
0066 bool Pythia6::Run() {
0067   TPythia6* pythia = TPythia6::Instance();
0068   TStopwatch timer;
0069   double lastTime(0.);
0070   while (mTree->GetEntries() < mNEvents) {
0071     const int initialNGenerated = pythia->GetMSTI(5);
0072     const int initialNTrials = pythia->GetPyint5()->NGEN[2][0];
0073     TBranch* branch = mTree->GetBranch("event");
0074     mFactory->Fill(*branch);
0075     mTree->Fill();
0076     // Count the number of generations and trials for this event
0077     mNGenerated += pythia->GetMSTI(5) - initialNGenerated;
0078     const int trials = pythia->GetPyint5()->NGEN[2][0] - initialNTrials;
0079     mNTrials += trials;
0080     if ((mTree->GetEntries() % mPrintInterval) == 0) {
0081       double time = timer.RealTime();
0082       std::cout << mTree->GetEntries() << " events in " <<
0083       time << " seconds (+" << time - lastTime << ")" << std::endl;
0084       lastTime = time;
0085       timer.Start(false);
0086     }  // if
0087   }  // while
0088   // Write the tree
0089   mFile->cd();
0090   mTree->Write();
0091   mFile->Purge();
0092   // Write run information
0093   std::stringstream ss;
0094   std::string s;
0095   // Write the total cross section
0096   ss << pythia->GetPARI(1) * 1000.;  // * 1000 --> microbarn
0097   ss >> s;
0098   TObjString(s.c_str()).Write("crossSection");
0099   // Write the number of generated events
0100   ss.str("");
0101   ss.clear();
0102   ss << mNGenerated;
0103   ss >> s;
0104   TObjString(s.c_str()).Write("nEvents");
0105   // Write the number of trials required to make the events
0106   ss.str("");
0107   ss.clear();
0108   ss << mNTrials;
0109   ss >> s;
0110   TObjString(s.c_str()).Write("nTrials");
0111   return true;
0112 }
0113
0114 }  // namespace erhic