Notification/src/TmpSimEvent.cc

0001 #include "SimG4Core/Notification/interface/TmpSimEvent.h"
0002 #include "SimDataFormats/EncodedEventId/interface/EncodedEventId.h"
0003
0004 #include <CLHEP/Units/SystemOfUnits.h>
0005
0006 class IdSort {
0007 public:
0008   bool operator()(const SimTrack& a, const SimTrack& b) { return a.trackId() < b.trackId(); }
0009 };
0010
0011 TmpSimEvent::TmpSimEvent() {
0012   g4vertices_.reserve(2000);
0013   g4tracks_.reserve(4000);
0014 }
0015
0016 TmpSimEvent::~TmpSimEvent() { clear(); }
0017
0018 void TmpSimEvent::clear() {
0019   for (auto& ptr : g4tracks_) {
0020     delete ptr;
0021   }
0022   g4tracks_.clear();
0023   for (auto& ptr : g4vertices_) {
0024     delete ptr;
0025   }
0026   g4vertices_.clear();
0027 }
0028
0029 void TmpSimEvent::load(edm::SimTrackContainer& c) const {
0030   const double invgev = 1.0 / CLHEP::GeV;
0031   for (auto& trk : g4tracks_) {
0032     int ip = trk->part();
0033     const math::XYZVectorD& mom = trk->momentum();
0034     math::XYZTLorentzVectorD p(mom.x() * invgev, mom.y() * invgev, mom.z() * invgev, trk->energy() * invgev);
0035     int iv = trk->ivert();
0036     int ig = trk->igenpart();
0037     int id = trk->id();
0038     bool isBackScatter = trk->isFromBackScattering();
0039     bool isPrimary = trk->isPrimary();
0040     int primaryGenPartId = trk->getPrimaryID();  // filled if the G4Track had this info
0041     // ip = particle ID as PDG
0042     // p  = 4-momentum in GeV
0043     // iv = corresponding TmpSimVertex index
0044     // ig = corresponding GenParticle index only if is a primary, otherwise is -1
0045     // primaryGenPartId = corresponding GenParticle index also if not primary
0046     // id = corresponding g4Track Id
0047     SimTrack t = SimTrack(ip, p, iv, ig, trk->trackerSurfacePosition(), trk->trackerSurfaceMomentum());
0048     t.setTrackId(id);
0049     t.setEventId(EncodedEventId(0));
0050     t.setCrossedBoundaryVars(
0051         trk->crossedBoundary(), trk->getIDAtBoundary(), trk->getPositionAtBoundary(), trk->getMomentumAtBoundary());
0052     if (isBackScatter)
0053       t.setFromBackScattering();
0054     if (isPrimary)
0055       t.setIsPrimary();
0056     else
0057       t.setGenParticleID(primaryGenPartId);
0058     c.push_back(t);
0059   }
0060   std::stable_sort(c.begin(), c.end(), IdSort());
0061 }
0062
0063 void TmpSimEvent::load(edm::SimVertexContainer& c) const {
0064   const double invcm = 1.0 / CLHEP::cm;
0065   // index of the vertex is needed to make SimVertex object
0066   for (unsigned int i = 0; i < g4vertices_.size(); ++i) {
0067     TmpSimVertex* vtx = g4vertices_[i];
0068     auto pos = vtx->vertexPosition();
0069     math::XYZVectorD v3(pos.x() * invcm, pos.y() * invcm, pos.z() * invcm);
0070     float t = (float)(vtx->vertexGlobalTime() / CLHEP::second);
0071     int iv = vtx->parentIndex();
0072     // v3 = position in cm
0073     // t  = global time in second
0074     // iv = index of the parent in the SimEvent SimTrack container (-1 if no parent)
0075     SimVertex v = SimVertex(v3, t, iv, i);
0076     v.setProcessType((unsigned int)vtx->processType());
0077     v.setEventId(EncodedEventId(0));
0078     c.push_back(v);
0079   }
0080 }