Notification/src/SimTrackManager.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:     Application
0004 // Class  :     SimTrackManager
0005 //
0006 // Implementation:
0007 //     <Notes on implementation>
0008 //
0009 // Original Author:
0010 //         Created:  Fri Nov 25 17:44:19 EST 2005
0011 //
0012
0013 // system include files
0014 #include <iostream>
0015
0016 // user include files
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 //#define DebugLog
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 /// this saves a track and all its parents looping over the non ordered vector
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   // fill the map with the final mother-daughter relationship
0106   idsave.swap(ancestorList);
0107   std::stable_sort(idsave.begin(), idsave.end());
0108
0109   std::vector<std::pair<int, int> >().swap(ancestorList);
0110
0111   // to get a backward compatible order
0112   std::stable_sort(m_trackContainer.begin(), m_trackContainer.end(), trkIDLess());
0113
0114   // to reset the GenParticle ID of a SimTrack to its pre-LHCTransport value
0115   resetGenID();
0116
0117   reallyStoreTracks(simEvent);
0118 }
0119
0120 void SimTrackManager::reallyStoreTracks(TmpSimEvent* simEvent) {
0121   // loop over the (now ordered) vector and really save the tracks
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     // at this stage there is one vertex per track,
0129     // so the vertex id of track N is also N
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     // Track at surface is the track at intersection point between tracker and calo
0144     // envelops if exist. If not exist the momentum is zero, position is the end of
0145     // the track
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     // loop over saved vertices
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       // first loop forward
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       // recursive loop
0227       do {
0228         notFound = true;
0229         // search ID scan backward
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           // ID not in the list of saved track - look into ancestors
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 }