File indexing completed on 2023-04-20 00:29:47
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 #include <iostream>
0015
0016
0017 #include "SimG4Core/Notification/interface/SimTrackManager.h"
0018 #include "SimG4Core/Notification/interface/TmpSimTrack.h"
0019 #include "SimG4Core/Notification/interface/TmpSimVertex.h"
0020 #include "SimG4Core/Notification/interface/TmpSimEvent.h"
0021 #include "SimG4Core/Notification/interface/TrackInformation.h"
0022
0023 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0024 #include "FWCore/Utilities/interface/Exception.h"
0025
0026 #include "G4VProcess.hh"
0027 #include "G4Track.hh"
0028 #include "G4ThreeVector.hh"
0029 #include "G4SystemOfUnits.hh"
0030
0031
0032
0033 SimTrackManager::SimTrackManager() {
0034 idsave.reserve(1000);
0035 ancestorList.reserve(1000);
0036 m_trackContainer.reserve(1000);
0037 m_endPoints.reserve(1000);
0038 }
0039
0040 SimTrackManager::~SimTrackManager() { reset(); }
0041
0042 void SimTrackManager::reset() {
0043 deleteTracks();
0044 cleanVertexMap();
0045 idsave.clear();
0046 ancestorList.clear();
0047 lastTrack = 0;
0048 lastHist = 0;
0049 }
0050
0051 void SimTrackManager::deleteTracks() {
0052 if (!m_trackContainer.empty()) {
0053 for (auto& ptr : m_trackContainer) {
0054 delete ptr;
0055 }
0056 m_trackContainer.clear();
0057 m_endPoints.clear();
0058 }
0059 }
0060
0061 void SimTrackManager::cleanVertexMap() {
0062 m_vertexMap.clear();
0063 m_vertexMap.swap(m_vertexMap);
0064 m_nVertices = 0;
0065 }
0066
0067 void SimTrackManager::addTrack(TrackWithHistory* iTrack, const G4Track* track, bool inHistory, bool withAncestor) {
0068 std::pair<int, int> thePair(iTrack->trackID(), iTrack->parentID());
0069 idsave.push_back(thePair);
0070 if (inHistory) {
0071 m_trackContainer.push_back(iTrack);
0072 const auto& v = track->GetStep()->GetPostStepPoint()->GetPosition();
0073 const double invcm = 1.0 / CLHEP::cm;
0074 std::pair<int, math::XYZVectorD> p(iTrack->trackID(),
0075 math::XYZVectorD(v.x() * invcm, v.y() * invcm, v.z() * invcm));
0076 m_endPoints.push_back(p);
0077 }
0078 if (withAncestor) {
0079 std::pair<int, int> thisPair(iTrack->trackID(), 0);
0080 ancestorList.push_back(thisPair);
0081 }
0082 }
0083
0084
0085 void SimTrackManager::saveTrackAndItsBranch(TrackWithHistory* trkWHist) {
0086 TrackWithHistory* trkH = trkWHist;
0087 if (trkH == nullptr) {
0088 edm::LogError("SimTrackManager") << " SimTrackManager::saveTrackAndItsBranch got 0 pointer ";
0089 throw cms::Exception("SimTrackManager::saveTrackAndItsBranch") << " cannot handle hits for tracking";
0090 }
0091 trkH->setToBeSaved();
0092 unsigned int parent = trkH->parentID();
0093
0094 auto tk_itr =
0095 std::lower_bound(m_trackContainer.begin(), m_trackContainer.end(), parent, SimTrackManager::StrictWeakOrdering());
0096
0097 if (tk_itr != m_trackContainer.end() && (*tk_itr)->trackID() == parent) {
0098 saveTrackAndItsBranch(*tk_itr);
0099 }
0100 }
0101
0102 void SimTrackManager::storeTracks(TmpSimEvent* simEvent) {
0103 cleanTracksWithHistory();
0104
0105
0106 idsave.swap(ancestorList);
0107 std::stable_sort(idsave.begin(), idsave.end());
0108
0109 std::vector<std::pair<int, int> >().swap(ancestorList);
0110
0111
0112 std::stable_sort(m_trackContainer.begin(), m_trackContainer.end(), trkIDLess());
0113
0114
0115 resetGenID();
0116
0117 reallyStoreTracks(simEvent);
0118 }
0119
0120 void SimTrackManager::reallyStoreTracks(TmpSimEvent* simEvent) {
0121
0122 #ifdef DebugLog
0123 edm::LogVerbatim("SimTrackManager") << "reallyStoreTracks() NtracksWithHistory= " << m_trackContainer->size();
0124 #endif
0125
0126 int nn = m_endPoints.size();
0127 for (auto& trkH : m_trackContainer) {
0128
0129
0130 unsigned int iParentID = trkH->parentID();
0131 int ig = trkH->genParticleID();
0132 int ivertex = getOrCreateVertex(trkH, iParentID, simEvent);
0133
0134 auto ptr = trkH;
0135 if (0 < iParentID) {
0136 for (auto& trk : m_trackContainer) {
0137 if (trk->trackID() == iParentID) {
0138 ptr = trk;
0139 break;
0140 }
0141 }
0142 }
0143
0144
0145
0146 const math::XYZVectorD& pm = ptr->momentum();
0147 math::XYZVectorD spos(0., 0., 0.);
0148 math::XYZTLorentzVectorD smom(0., 0., 0., 0.);
0149 int id = trkH->trackID();
0150 if (trkH->crossedBoundary()) {
0151 spos = trkH->getPositionAtBoundary();
0152 smom = trkH->getMomentumAtBoundary();
0153 } else {
0154 for (int i = 0; i < nn; ++i) {
0155 if (id == m_endPoints[i].first) {
0156 spos = m_endPoints[i].second;
0157 break;
0158 }
0159 }
0160 }
0161
0162 TmpSimTrack* g4simtrack =
0163 new TmpSimTrack(id, trkH->particleID(), trkH->momentum(), trkH->totalEnergy(), ivertex, ig, pm, spos, smom);
0164 g4simtrack->copyCrossedBoundaryVars(trkH);
0165 simEvent->add(g4simtrack);
0166 }
0167 }
0168
0169 int SimTrackManager::getOrCreateVertex(TrackWithHistory* trkH, int iParentID, TmpSimEvent* simEvent) {
0170 int parent = -1;
0171 for (auto& trk : m_trackContainer) {
0172 int id = trk->trackID();
0173 if (id == iParentID) {
0174 parent = id;
0175 break;
0176 }
0177 }
0178
0179 VertexMap::const_iterator iterator = m_vertexMap.find(parent);
0180 if (iterator != m_vertexMap.end()) {
0181
0182 for (auto& xx : m_vertexMap[parent]) {
0183 if ((trkH->vertexPosition() - xx.second).Mag2() < 1.e-6) {
0184 return xx.first;
0185 }
0186 }
0187 }
0188
0189 unsigned int ptype = 0;
0190 const G4VProcess* pr = trkH->creatorProcess();
0191 if (nullptr != pr) {
0192 ptype = pr->GetProcessSubType();
0193 }
0194 simEvent->add(new TmpSimVertex(trkH->vertexPosition(), trkH->globalTime(), parent, ptype));
0195 m_vertexMap[parent].push_back(VertexPosition(m_nVertices, trkH->vertexPosition()));
0196 ++m_nVertices;
0197 return (m_nVertices - 1);
0198 }
0199
0200 int SimTrackManager::idSavedTrack(int id) const {
0201 int idMother = id;
0202 if (id > 0) {
0203 unsigned int n = idsave.size();
0204 if (0 < n) {
0205 int jmax = n - 1;
0206 int j, id1;
0207
0208
0209 bool notFound = true;
0210 for (j = 0; j <= jmax; ++j) {
0211 if ((idsave[j]).first == idMother) {
0212 id1 = (idsave[j]).second;
0213 if (0 == id1 || id1 == idMother) {
0214 return id1;
0215 }
0216 jmax = j - 1;
0217 idMother = id1;
0218 notFound = false;
0219 break;
0220 }
0221 }
0222 if (notFound) {
0223 return 0;
0224 }
0225
0226
0227 do {
0228 notFound = true;
0229
0230 for (j = jmax; j >= 0; --j) {
0231 if ((idsave[j]).first == idMother) {
0232 id1 = (idsave[j]).second;
0233 if (0 == id1 || id1 == idMother) {
0234 return id1;
0235 }
0236 jmax = j - 1;
0237 idMother = id1;
0238 notFound = false;
0239 break;
0240 }
0241 }
0242 if (notFound) {
0243
0244 jmax = ancestorList.size() - 1;
0245 for (j = jmax; j >= 0; --j) {
0246 if ((ancestorList[j]).first == idMother) {
0247 idMother = (ancestorList[j]).second;
0248 return idMother;
0249 }
0250 }
0251 return 0;
0252 }
0253 } while (!notFound);
0254 }
0255 }
0256 return idMother;
0257 }
0258
0259 void SimTrackManager::fillMotherList() {
0260 if (!ancestorList.empty() && lastHist > ancestorList.size()) {
0261 lastHist = ancestorList.size();
0262 edm::LogError("SimTrackManager") << " SimTrackManager::fillMotherList track index corrupted";
0263 }
0264 #ifdef DebugLog
0265 edm::LogVerbatim("SimTrackManager") << "### SimTrackManager::fillMotherList: " << idsave.size()
0266 << " saved; ancestor: " << lastHist << " " << ancestorList.size();
0267 for (unsigned int i = 0; i < idsave.size(); ++i) {
0268 edm::LogVerbatim("SimTrackManager") << " ISV: Track ID = " << (idsave[i]).first
0269 << " Mother ID = " << (idsave[i]).second;
0270 }
0271 #endif
0272 for (unsigned int n = lastHist; n < ancestorList.size(); ++n) {
0273 int theMotherId = idSavedTrack((ancestorList[n]).first);
0274 ancestorList[n].second = theMotherId;
0275 #ifdef DebugLog
0276 LogDebug("SimTrackManager") << "Track ID = " << (ancestorList[n]).first
0277 << " Mother ID = " << (ancestorList[n]).second;
0278 #endif
0279 }
0280
0281 lastHist = ancestorList.size();
0282 idsave.clear();
0283 }
0284
0285 void SimTrackManager::cleanTracksWithHistory() {
0286 if (m_trackContainer.empty() && idsave.empty()) {
0287 return;
0288 }
0289
0290 #ifdef DebugLog
0291 LogDebug("SimTrackManager") << "SimTrackManager::cleanTracksWithHistory has " << idsave.size()
0292 << " mother-daughter relationships stored with lastTrack = " << lastTrack;
0293 #endif
0294
0295 if (lastTrack > 0 && lastTrack >= m_trackContainer.size()) {
0296 lastTrack = 0;
0297 edm::LogError("SimTrackManager") << " SimTrackManager::cleanTracksWithHistory track index corrupted";
0298 }
0299
0300 std::stable_sort(m_trackContainer.begin() + lastTrack, m_trackContainer.end(), trkIDLess());
0301 std::stable_sort(idsave.begin(), idsave.end());
0302
0303 #ifdef DebugLog
0304 LogDebug("SimTrackManager") << " SimTrackManager::cleanTracksWithHistory knows " << m_trksForThisEvent->size()
0305 << " tracks with history before branching";
0306 for (unsigned int it = 0; it < m_trackContainer.size(); it++) {
0307 LogDebug("SimTrackManager") << " 1 - Track in position " << it << " G4 track number "
0308 << m_trackContainer[it]->trackID() << " mother " << m_trackContainer[it]->parentID()
0309 << " status " << m_trackContainer[it]->saved();
0310 }
0311 #endif
0312
0313 for (auto& t : m_trackContainer) {
0314 if (t->saved()) {
0315 saveTrackAndItsBranch(t);
0316 }
0317 }
0318 unsigned int num = lastTrack;
0319 for (unsigned int it = lastTrack; it < m_trackContainer.size(); ++it) {
0320 auto t = m_trackContainer[it];
0321 int g4ID = m_trackContainer[it]->trackID();
0322 if (t->saved()) {
0323 if (it > num) {
0324 m_trackContainer[num] = t;
0325 }
0326 ++num;
0327 for (auto& xx : idsave) {
0328 if (xx.first == g4ID) {
0329 xx.second = g4ID;
0330 break;
0331 }
0332 }
0333 } else {
0334 delete t;
0335 }
0336 }
0337
0338 m_trackContainer.resize(num);
0339
0340 #ifdef DebugLog
0341 LogDebug("SimTrackManager") << " AFTER CLEANING, I GET " << m_trackContainer.size()
0342 << " tracks to be saved persistently";
0343 for (unsigned int it = 0; it < m_trackContainer.size(); ++it) {
0344 LogDebug("SimTrackManager") << " Track in position " << it << " G4 track number " << m_trackContainer[it]->trackID()
0345 << " mother " << m_trackContainer[it]->parentID() << " Status "
0346 << m_trackContainer[it]->saved() << " id " << m_trackContainer[it]->particleID()
0347 << " E(MeV)= " << m_trackContainer[it]->totalEnergy();
0348 }
0349 #endif
0350
0351 fillMotherList();
0352 lastTrack = m_trackContainer.size();
0353 }
0354
0355 void SimTrackManager::resetGenID() {
0356 if (theLHCTlink == nullptr)
0357 return;
0358
0359 for (auto& trkH : m_trackContainer) {
0360 int genParticleID = trkH->genParticleID();
0361 if (genParticleID == -1) {
0362 continue;
0363 } else {
0364 for (auto& xx : *theLHCTlink) {
0365 if (xx.afterHector() == genParticleID) {
0366 trkH->setGenParticleID(xx.beforeHector());
0367 continue;
0368 }
0369 }
0370 }
0371 }
0372
0373 theLHCTlink = nullptr;
0374 }
0375
0376 void SimTrackManager::ReportException(unsigned int id) const {
0377 throw cms::Exception("Unknown", "SimTrackManager::getTrackByID")
0378 << "Fail to get track " << id << " from SimTrackManager, container size= " << m_trackContainer.size();
0379 }