src/erhic/EventMC.cxx

0001 /**
0002  \file
0003  Implementation of class erhic::EventMC.
0004
0005  \author    Thomas Burton
0006  \date      2011-10-31
0007  \copyright 2011 Brookhaven National Lab
0008  */
0009
0010 #include "eicsmear/erhic/EventMC.h"
0011
0012 #include <iostream>
0013 #include <list>
0014 #include <string>
0015 #include <vector>
0016
0017 #include <TCollection.h>
0018 #include <TDatabasePDG.h>
0019 #include <TDirectory.h>
0020 #include <TParticlePDG.h>
0021 #include <TTree.h>
0022
0023 #include "eicsmear/erhic/ParticleIdentifier.h"
0024 #include "eicsmear/erhic/ParticleMC.h"
0025
0026 namespace erhic {
0027
0028 // We use these vectors/iterators a lot, so define
0029 // some typedefs for convenience.
0030 typedef std::vector<const VirtualParticle*> TrackVector;
0031 typedef TrackVector::iterator TrackVectorIter;
0032 typedef TrackVector::const_iterator TrackVectorCIter;
0033
0034 EventMC::EventMC()
0035 : number(-1)
0036 , process(-1)
0037 , nTracks(-1)
0038 , ELeptonInNucl(NAN)
0039 , ELeptonOutNucl(NAN)
0040 , particles("erhic::ParticleMC", 100) {
0041 }
0042
0043 EventMC::~EventMC() {
0044   // No memory to clear. The TClonesArray takes care of the
0045   // particles when it is destroyed.
0046 }
0047
0048 TrackVector EventMC::GetTracks() const {
0049   TrackVector tracks;
0050   TObject* object(NULL);
0051   TIter next(&particles);
0052   while ((object = next())) {
0053     tracks.push_back(static_cast<ParticleMC*>(object));
0054   }  // while
0055   return tracks;
0056 }
0057
0058 void EventMC::HadronicFinalState(TrackVector& final_) const {
0059   // This is simple but not very efficient.
0060   // Copy vector to a list and use list::remove to get rid
0061   // of the scattered lepton. Then copy this list back to
0062   // the vector.
0063   // Note that the method is a bit of a misnomer - it will return ALL final
0064   // particles other than the scattered lepton
0065   // (intentionally, since you want to take decay products into account as well)
0066   FinalState(final_);
0067   std::list<const VirtualParticle*> plist(final_.begin(),
0068                                           final_.end());
0069   plist.remove(ScatteredLepton());
0070   final_ = TrackVector(plist.begin(), plist.end());
0071 }
0072
0073 // Get the particles that belong to the hadronic final state.
0074 // The stored Particle* are pointers to the original particles in the event
0075 // so don't delete them!
0076 void EventMC::FinalState(TrackVector& final_) const {
0077   final_.clear();
0078   TIter next(&particles);
0079   ParticleMC* p(NULL);
0080   while ((p = static_cast<ParticleMC*>(next()))) {
0081     if (1 == p->GetStatus()) {
0082       final_.push_back(p);
0083     }  // if
0084   }  // while
0085 }
0086
0087 TLorentzVector EventMC::FinalStateMomentum() const {
0088   TrackVector final_;
0089   FinalState(final_);
0090   TLorentzVector mom;
0091   for (TrackVectorCIter i = final_.begin(); i != final_.end(); ++i) {
0092     mom += (*i)->Get4Vector();
0093   }  // for
0094
0095   return mom;
0096 }
0097
0098 TLorentzVector EventMC::HadronicFinalStateMomentum() const {
0099   TrackVector final_;
0100   HadronicFinalState(final_);
0101   TLorentzVector mom;
0102   for (TrackVectorCIter i = final_.begin(); i != final_.end(); ++i) {
0103     mom += (*i)->Get4Vector();
0104   }  // for
0105
0106   return mom;
0107 }
0108
0109 Double_t EventMC::FinalStateCharge() const {
0110   TrackVector final_;
0111   FinalState(final_);
0112   Double_t charge(0);
0113   TDatabasePDG* pdg = TDatabasePDG::Instance();
0114   for (TrackVectorCIter i = final_.begin(); i != final_.end(); ++i) {
0115     TParticlePDG* part = pdg->GetParticle((*i)->Id());
0116     if (part) {
0117       charge += part->Charge() / 3.;
0118     } else {
0119       std::cout << "Unknown particle: " << (*i)->Id() << std::endl;
0120     }  // if
0121   }  // for
0122   return charge;
0123 }
0124
0125   // See header!
0126 const ParticleMC* EventMC::BeamLepton() const {
0127   return GetTrack(0);
0128 }
0129
0130   // See header!
0131 const ParticleMC* EventMC::BeamHadron() const {
0132   return GetTrack(1);
0133 }
0134
0135   // See header!
0136 const ParticleMC* EventMC::ScatteredLepton() const {
0137   return GetTrack(2);
0138 }
0139
0140   // See header!
0141 const ParticleMC* EventMC::ExchangeBoson() const {
0142   return GetTrack(3);
0143 }
0144
0145
0146 void EventMC::Reset() {
0147   number = -1;
0148   process = -1;
0149   nTracks = -1;
0150   x = QSquared = y = WSquared = nu = ELeptonInNucl = ELeptonOutNucl = NAN;
0151 }
0152
0153 void EventMC::Clear(Option_t* /*option*/) {
0154   Reset();
0155   particles.Clear();
0156 }
0157
0158 void EventMC::AddLast(ParticleMC* track) {
0159   new(particles[particles.GetEntries()]) ParticleMC(*track);
0160   nTracks = particles.GetEntries();
0161 }
0162
0163 void EventMC::Print( const Option_t *option) const {
0164   EventDis::Print();
0165   std::cout << "I \t KS \t id \t orig\t daughter \t ldaughter \t "
0166     << " px \t py \t pz \t E \t m \t xv \t yv \t zv"
0167     << std::endl;
0168
0169        for( unsigned int t=0; t<GetNTracks(); ++t) {
0170      const Particle* inParticle = GetTrack(t);
0171      inParticle->Print();
0172        }
0173 }
0174
0175   //
0176 // class Reader
0177 //
0178
0179 Reader::Reader(const std::string& treeName)
0180 : mEvent(NULL)
0181 , mTree(NULL) {
0182   if (gDirectory) gDirectory->GetObject(treeName.c_str(), mTree);
0183   if (mTree) mTree->SetBranchAddress("event", &mEvent);
0184 }
0185
0186 EventMC* Reader::Read(Long64_t i) {
0187   EventMC* event(NULL);
0188   if (mTree) {
0189     mTree->GetEntry(i);
0190     event = mEvent;
0191   }  // if
0192   return event;
0193 }
0194
0195 }  // namespace erhic