src/erhic/EventFactory.cxx

0001 /**
0002    \file
0003    Implementation of class erhic::EventFactory.
0004
0005    \author    Thomas Burton
0006    \date      2011-10-31
0007    \copyright 2011 Brookhaven National Lab
0008 */
0009
0010 #include "eicsmear/erhic/EventFactory.h"
0011
0012 #include <memory>
0013 #include <stdexcept>
0014 #include <string>
0015
0016 #include <TClass.h>
0017 #include <TProcessID.h>
0018
0019 #include "eicsmear/erhic/BeamParticles.h"
0020 #include "eicsmear/erhic/EventPythia.h"
0021 #include "eicsmear/erhic/EventHepMC.h"
0022 #include "eicsmear/erhic/EventMilou.h"
0023 #include "eicsmear/erhic/EventDjangoh.h"
0024 #include "eicsmear/erhic/EventDpmjet.h"
0025 #include "eicsmear/erhic/EventRapgap.h"
0026 #include "eicsmear/erhic/EventPepsi.h"
0027 #include "eicsmear/erhic/EventGmcTrans.h"
0028 #include "eicsmear/erhic/EventSimple.h"
0029 #include "eicsmear/erhic/EventDEMP.h"
0030 #include "eicsmear/erhic/EventSartre.h"
0031 #include "eicsmear/functions.h"  // For getFirstNonBlank()
0032 #include "eicsmear/erhic/Kinematics.h"
0033 #include "eicsmear/erhic/ParticleIdentifier.h"
0034 #include "eicsmear/erhic/ParticleMC.h"
0035
0036 #include <TVector3.h>
0037 #include <TParticlePDG.h>
0038 #include <TLorentzVector.h>
0039 #include <TDatabasePDG.h>
0040
0041 #include<map>
0042
0043 using std::cout;
0044 using std::cerr;
0045 using std::endl;
0046 using std::map;
0047
0048 namespace erhic {
0049
0050   template<typename T>
0051   bool EventFromAsciiFactory<T>::AtEndOfEvent() const {
0052     return !mLine.empty()
0053       && mLine.find("finished") != std::string::npos;
0054   }
0055
0056   // Use this struct to automatically reset TProcessID object count.
0057   struct TProcessIdObjectCount {
0058     // Initialse object with current TProcessID object count.
0059     TProcessIdObjectCount() {
0060       count = TProcessID::GetObjectCount();
0061     }
0062     // Restore object count to the value at initialisation.
0063     // See example in $ROOTSYS/test/Event.cxx
0064     // To save space in the table keeping track of all referenced objects
0065     // we assume that our events do not address each other.
0066     ~TProcessIdObjectCount() {
0067       TProcessID::SetObjectCount(count);
0068     }
0069     int count;
0070   };
0071
0072   template<typename T>
0073   T* EventFromAsciiFactory<T>::Create() {
0074     // Save current object count. Will reset it when this function returns.
0075     TProcessIdObjectCount objectCount;
0076     mEvent.reset(new T);
0077     // We use this flag to check input doesn't end mid-event.
0078     // Initialised finished flag to "success" in case of no input.
0079     int finished(0);
0080     std::string error;
0081     // Read line-by-line until the stream is not good or we break out.
0082     while (std::getline(*mInput, mLine).good()) {
0083       // Reached end-of-event marker
0084       if (AtEndOfEvent()) {
0085     // If we built a good event (i.e. no errors reading strings)
0086     // perform end-of-event calculations.
0087     if (!error.empty()) {
0088       throw std::runtime_error(error);
0089     } else {
0090       finished = FinishEvent();  // 0 upon success
0091       break;
0092     }  // if
0093       } else if ('0' == getFirstNonBlank(mLine)) {
0094     // '0' indicates the event header line
0095     // An event started, set finished flag to "unfinished".
0096     finished = -1;
0097     // Parse string and check for validity.
0098     if (!mEvent->Parse(mLine)) {
0099       // Set error message based on bad event input.
0100       // Don't break out of the loop yet, so that we continue reading
0101       // lines until the end-of-event marker. That way we stay
0102       // "aligned" with the input data ready for the next event.
0103       error = "Bad event input: " + mLine;
0104     }  // if
0105       } else if ('=' != getFirstNonBlank(mLine)) {
0106     // Anything remaining other than a line of '=' is a particle line
0107     // Don't raise an exception for a failed track, as the event in
0108     // general may be OK. AddParticle will print a message though.
0109     if (!AddParticle()) {
0110       error = "Bad particle input in event";
0111     }  // if
0112       }  // if
0113     }  // if
0114     // Check for events that started but did not finish.
0115     // We should have ended with an end-of-event marker, meaning the finished
0116     // flag is set to zero. If not, the finished flag will be non-zero.
0117     if (finished != 0) {
0118       mEvent.reset(NULL);
0119       throw std::runtime_error("Ended mid-event");
0120     }  // if
0121     // Return a NULL event *without* throwing an exception to indicate
0122     // end-of-file. We shouldn't have hit eof if we have read a good event,
0123     // as we won't have yet tried to read past the end (mLine will still be
0124     // at the end-of-file marker line).
0125     if (mInput->eof()) {
0126       mEvent.reset(NULL);
0127     }  // if
0128     return mEvent.release();
0129   }
0130
0131   template<typename T>
0132   Int_t EventFromAsciiFactory<T>::FinishEvent() {
0133     // First, find the beams, exchange boson, scattered lepton.
0134     // if this isn't a DIS event, skip Kinematics computation
0135     BeamParticles beams;
0136     if (!ParticleIdentifier::IdentifyBeams(*mEvent, beams)) {
0137       return 0;  // must not signal error.
0138       // std::cerr << "EventFromAsciiFactory::FinishEvent(): failed to find beams"  << std::endl;
0139       // return -1;
0140     }  // if
0141
0142     std::unique_ptr<DisKinematics> nm( LeptonKinematicsComputer(*mEvent).Calculate());
0143     std::unique_ptr<DisKinematics> jb( JacquetBlondelComputer(*mEvent).Calculate());
0144     std::unique_ptr<DisKinematics> da( DoubleAngleComputer(*mEvent).Calculate());
0145     if (nm.get()) {
0146       mEvent->SetLeptonKinematics(*nm);
0147     }  // if
0148     for (unsigned n(0); n < mEvent->GetNTracks(); ++n) {
0149       mEvent->GetTrack(n)->ComputeEventDependentQuantities(*mEvent);
0150     }  // for
0151     if (jb.get()) {
0152       mEvent->SetJacquetBlondelKinematics(*jb);
0153     }  // if
0154     if (da.get()) {
0155       mEvent->SetDoubleAngleKinematics(*da);
0156     }  // if
0157
0158     // We also have to set the remaining variables not taken care of
0159     // by the general DIS event kinematic computations.
0160     const TLorentzVector h = beams.BeamHadron();
0161     TLorentzVector l = beams.BeamLepton();
0162     TLorentzVector s = beams.ScatteredLepton();
0163     TVector3 boost = -h.BoostVector();
0164     l.Boost(boost);
0165     s.Boost(boost);
0166     mEvent->SetELeptonInNuclearFrame(l.E());
0167     mEvent->SetEScatteredInNuclearFrame(s.E());
0168     return 0;
0169   }
0170
0171   template<typename T>
0172   bool EventFromAsciiFactory<T>::AddParticle() {
0173     //return true;
0174     try {
0175       if (mEvent.get()) {
0176     ParticleMC particle(mLine, mEvent->RequiresEaParticleFields());  // Throws if the string is bad
0177     particle.SetEvent(mEvent.get());
0178     mEvent->AddLast(&particle);
0179     //ParticleMCeA *particle = new ParticleMCeA(mLine);  // Throws if the string is bad
0180     //particle->SetEvent(mEvent.get());
0181     //mEvent->AddLast(particle);
0182     //delete particle;
0183       }  // if
0184       return true;
0185     }  // try
0186     catch(std::exception& error) {
0187       std::cerr << "Exception building particle: " << error.what() << std::endl;
0188       return false;
0189     }
0190   }
0191
0192   template<typename T>
0193   std::string EventFromAsciiFactory<T>::EventName() const {
0194     return T::Class()->GetName();
0195   }
0196
0197   template<typename T>
0198   void EventFromAsciiFactory<T>::FindFirstEvent()  {
0199     for (int i=0; i<5; i++)
0200       {
0201     std::getline(*mInput,mLine);
0202       }
0203   }
0204
0205
0206 }  // namespace erhic
0207
0208 namespace {
0209
0210   // Need this to generate the CINT code for each version
0211   erhic::EventFromAsciiFactory<erhic::EventDjangoh> ed;
0212   erhic::EventFromAsciiFactory<erhic::EventDpmjet> ej;
0213   erhic::EventFromAsciiFactory<erhic::EventPepsi> ee;
0214   erhic::EventFromAsciiFactory<erhic::EventMilou> em;
0215   erhic::EventFromAsciiFactory<erhic::EventHepMC> eh;
0216   erhic::EventFromAsciiFactory<erhic::EventRapgap> er;
0217   erhic::EventFromAsciiFactory<erhic::EventPythia> ep;
0218   erhic::EventFromAsciiFactory<erhic::EventGmcTrans> eg;
0219   erhic::EventFromAsciiFactory<erhic::EventSimple> es;
0220   erhic::EventFromAsciiFactory<erhic::EventDEMP> edp;
0221   erhic::EventFromAsciiFactory<erhic::EventSartre> esa;
0222   erhic::EventFromAsciiFactory<erhic::EventBeagle> eb;
0223
0224 }  // namespace