src/smear/SmearTree.cxx

0001 /**
0002  \file
0003  Here lies the function for smearing trees
0004
0005  \author    Michael Savastio
0006  \date      2011-08-12
0007  \copyright 2011 Brookhaven National Lab
0008  */
0009
0010 #include <iostream>
0011 #include <fstream>
0012 #include <sstream>
0013 #include <cmath>
0014 #include <memory>
0015
0016 #include <TClass.h>
0017 #include <TROOT.h>
0018 #include <TSystem.h>
0019 #include <TChain.h>
0020 #include <TString.h>
0021 #include <TRandom2.h>
0022 #include <TTree.h>
0023 #include <TFile.h>
0024 #include <TStopwatch.h>
0025 #include <TH1D.h>
0026
0027 #include "eicsmear/erhic/VirtualParticle.h"
0028 #include "eicsmear/smear/Detector.h"
0029 #include "eicsmear/smear/EventDisFactory.h"
0030 #include "eicsmear/smear/ParticleMCS.h"
0031 #include "eicsmear/smear/Smear.h"
0032
0033 #ifdef WITH_PYTHIA6
0034 #include "eicsmear/hadronic/EventSmear.h"
0035 #endif
0036
0037 /**
0038  Smear nEvents events from the TTree named EICTree in the named input file,
0039  using the smearing definitions in the Detector.
0040  Write the resulting Smeared TTree to a file named outFileName.
0041  If nEvents <= 0 smear all events in the tree.
0042  Returns 0 upon success, 1 upon failure.
0043  */
0044 int SmearTree(const Smear::Detector& detector, const TString& inFileName,
0045               const TString& outFileName, Long64_t nEvents) {
0046   // Open the input file and get the Monte Carlo tree from it.
0047   // Complain and quit if we don't find the file or the tree.
0048   TFile inFile(inFileName, "READ");
0049   if (!inFile.IsOpen()) {
0050     std::cerr << "Unable to open " << inFileName << std::endl;
0051   }  // if
0052   TTree* mcTree(NULL);
0053   inFile.GetObject("EICTree", mcTree);
0054   if (!mcTree) {
0055     std::cerr << "Unable to find EICTree in " << inFileName << std::endl;
0056     return 1;
0057   }  // if
0058   std::unique_ptr<erhic::VirtualEventFactory> builder;
0059   // Need to determine the type of object in the tree to choose
0060   // the correct smeared event builder.
0061   TClass* branchClass = TClass::GetClass(mcTree->GetBranch("event")->GetClassName());
0062   if (branchClass->InheritsFrom("erhic::EventDis")) {
0063     builder.reset(new Smear::EventDisFactory(detector,
0064                                              *(mcTree->GetBranch("event"))));
0065 #ifdef WITH_PYTHIA6
0066   } else if (branchClass->InheritsFrom("erhic::hadronic::EventMC")) {
0067     builder.reset(new Smear::HadronicEventBuilder(detector,
0068                                              *(mcTree->GetBranch("event"))));
0069 #endif
0070   } else {
0071     std::cerr << branchClass->GetName() << " is not supported for smearing" <<
0072     std::endl;
0073   }  // if
0074   // Open the output file.
0075   // Complain and quit if something goes wrong.
0076   TString outName(outFileName);
0077   if (outName.IsNull()) {
0078     outName = TString(inFileName).ReplaceAll(".root", ".smear.root");
0079   }  // if
0080   TFile outFile(outName, "RECREATE");
0081   if (!outFile.IsOpen()) {
0082     std::cerr << "Unable to create " << outName << std::endl;
0083     return 1;
0084   }  // if
0085   TTree smearedTree("Smeared", "A tree of smeared Monte Carlo events");
0086   TBranch* eventbranch = builder->Branch(smearedTree, "eventS");
0087   if (mcTree->GetEntries() < nEvents || nEvents < 1) {
0088     nEvents = mcTree->GetEntries();
0089   }  // if
0090   std::cout <<
0091   "/-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-/"
0092   << std::endl;
0093   std::cout <<
0094   "/  Commencing Smearing of " << nEvents << " events."
0095   << std::endl;
0096   std::cout <<
0097   "/-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-/"
0098   << std::endl;
0099   for (Long64_t i(0); i < nEvents; i++) {
0100     if (i % 10000 == 0 && i != 0) {
0101       std::cout << "Processing event " << i << std::endl;
0102     }  // if
0103     mcTree->GetEntry(i);
0104     builder->Fill(*eventbranch);
0105   }  // for
0106   smearedTree.Write();
0107   detector.Write("detector");
0108   outFile.Purge();
0109   std::cout <<
0110   "|~~~~~~~~~~~~~~~~~~ Completed Successfully ~~~~~~~~~~~~~~~~~~~|"
0111   << std::endl;
0112   return 0;
0113 }