Notification/src/G4SimEvent.cc

0001 #include "SimG4Core/Notification/interface/G4SimEvent.h"
0002 #include "SimDataFormats/EncodedEventId/interface/EncodedEventId.h"
0003
0004 #include "G4SystemOfUnits.hh"
0005
0006 class IdSort {
0007 public:
0008   bool operator()(const SimTrack& a, const SimTrack& b) { return a.trackId() < b.trackId(); }
0009 };
0010
0011 G4SimEvent::G4SimEvent()
0012     : hepMCEvent_(nullptr),
0013       weight_(0),
0014       collisionPoint_(math::XYZTLorentzVectorD(0., 0., 0., 0.)),
0015       nparam_(0),
0016       param_(0) {
0017   g4vertices_.reserve(2000);
0018   g4tracks_.reserve(4000);
0019 }
0020
0021 G4SimEvent::~G4SimEvent() { clear(); }
0022
0023 void G4SimEvent::clear() {
0024   // per suggestion by Chris Jones, it's faster
0025   // that delete back() and pop_back()
0026   for (auto& ptr : g4tracks_) {
0027     delete ptr;
0028   }
0029   g4tracks_.clear();
0030   for (auto& ptr : g4vertices_) {
0031     delete ptr;
0032   }
0033   g4vertices_.clear();
0034 }
0035
0036 void G4SimEvent::load(edm::SimTrackContainer& c) const {
0037   for (auto& trk : g4tracks_) {
0038     int ip = trk->part();
0039     math::XYZTLorentzVectorD p(
0040         trk->momentum().x() / GeV, trk->momentum().y() / GeV, trk->momentum().z() / GeV, trk->energy() / GeV);
0041     int iv = trk->ivert();
0042     int ig = trk->igenpart();
0043     int id = trk->id();
0044     math::XYZVectorD tkpos(trk->trackerSurfacePosition().x() / cm,
0045                            trk->trackerSurfacePosition().y() / cm,
0046                            trk->trackerSurfacePosition().z() / cm);
0047     math::XYZTLorentzVectorD tkmom(trk->trackerSurfaceMomentum().x() / GeV,
0048                                    trk->trackerSurfaceMomentum().y() / GeV,
0049                                    trk->trackerSurfaceMomentum().z() / GeV,
0050                                    trk->trackerSurfaceMomentum().e() / GeV);
0051     // ip = particle ID as PDG
0052     // pp = 4-momentum
0053     // iv = corresponding G4SimVertex index
0054     // ig = corresponding GenParticle index
0055     SimTrack t = SimTrack(ip, p, iv, ig, tkpos, tkmom);
0056     t.setTrackId(id);
0057     t.setEventId(EncodedEventId(0));
0058     if (trk->crossedBoundary())
0059       t.setCrossedBoundaryVars(
0060           trk->crossedBoundary(), trk->getIDAtBoundary(), trk->getPositionAtBoundary(), trk->getMomentumAtBoundary());
0061     c.push_back(t);
0062   }
0063   std::stable_sort(c.begin(), c.end(), IdSort());
0064 }
0065
0066 void G4SimEvent::load(edm::SimVertexContainer& c) const {
0067   for (unsigned int i = 0; i < g4vertices_.size(); ++i) {
0068     G4SimVertex* vtx = g4vertices_[i];
0069     //
0070     // starting 1_1_0_pre3, SimVertex stores in cm !!!
0071     //
0072     math::XYZVectorD v3(vtx->vertexPosition().x() / cm, vtx->vertexPosition().y() / cm, vtx->vertexPosition().z() / cm);
0073     float t = vtx->vertexGlobalTime() / second;
0074     int iv = vtx->parentIndex();
0075     // vv = position
0076     // t  = global time
0077     // iv = index of the parent in the SimEvent SimTrack container (-1 if no parent)
0078     SimVertex v = SimVertex(v3, t, iv, i);
0079     v.setProcessType(vtx->processType());
0080     v.setEventId(EncodedEventId(0));
0081     c.push_back(v);
0082   }
0083 }