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