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 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0023
0024 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0025 #include "FWCore/Utilities/interface/Exception.h"
0026
0027 #include "G4VProcess.hh"
0028 #include "G4Track.hh"
0029 #include "G4ThreeVector.hh"
0030 #include <CLHEP/Units/SystemOfUnits.h>
0031
0032 //#define DebugLog
0033 namespace {
0034   const double invcm = 1.0 / CLHEP::cm;
0035   const double r_limit2 = 1.e-6;  // 10 micron in CMS units
0036 }  // namespace
0037
0038 SimTrackManager::SimTrackManager(TmpSimEvent* ptr, int) : m_simEvent(ptr) {
0039   idsave.reserve(1000);
0040   ancestorList.reserve(1000);
0041   m_trackContainer.reserve(1000);
0042   m_endPoints.reserve(1000);
0043 }
0044
0045 SimTrackManager::~SimTrackManager() { reset(); }
0046
0047 void SimTrackManager::reset() {
0048   deleteTracks();
0049   cleanVertexMap();
0050   idsave.clear();
0051   ancestorList.clear();
0052   lastTrack = 0;
0053   lastHist = 0;
0054 }
0055
0056 void SimTrackManager::deleteTracks() {
0057   if (!m_trackContainer.empty()) {
0058     for (auto const& ptr : m_trackContainer) {
0059       delete ptr;
0060     }
0061     m_trackContainer.clear();
0062     m_endPoints.clear();
0063   }
0064 }
0065
0066 void SimTrackManager::cleanVertexMap() {
0067   m_vertexMap.clear();
0068   m_vertexMap.swap(m_vertexMap);
0069   m_nVertices = 0;
0070 }
0071
0072 void SimTrackManager::addTrack(TrackWithHistory* iTrack, const G4Track* track, bool inHistory, bool withAncestor) {
0073   std::pair<int, int> thePair(iTrack->trackID(), iTrack->parentID());
0074   idsave.push_back(thePair);
0075   if (inHistory) {
0076     auto info = static_cast<const TrackInformation*>(track->GetUserInformation());
0077     if (info->isInTrkFromBackscattering())
0078       iTrack->setFromBackScattering();
0079     // set there for the *non-primary* tracks the genParticle ID associated with the G4Track
0080     // for the primaries this is done in the TrackWithHistory constructor.
0081     // In the constructor of TrackWithHistory the info isPrimary is saved and used
0082     // to give -1 if the track is not a primary.
0083     if (not iTrack->isPrimary())
0084       iTrack->setGenParticleID(info->mcTruthID());
0085     m_trackContainer.push_back(iTrack);
0086     const auto& v = track->GetStep()->GetPostStepPoint()->GetPosition();
0087     std::pair<int, math::XYZVectorD> p(iTrack->trackID(),
0088                                        math::XYZVectorD(v.x() * invcm, v.y() * invcm, v.z() * invcm));
0089     m_endPoints.push_back(p);
0090   }
0091   if (withAncestor) {
0092     std::pair<int, int> thisPair(iTrack->trackID(), 0);
0093     ancestorList.push_back(thisPair);
0094   }
0095 }
0096
0097 /// this saves a track and all its parents looping over the non ordered vector
0098 void SimTrackManager::saveTrackAndItsBranch(TrackWithHistory* trkWHist) {
0099   TrackWithHistory* trkH = trkWHist;
0100   if (trkH == nullptr) {
0101     edm::LogError("SimTrackManager") << " SimTrackManager::saveTrackAndItsBranch got 0 pointer ";
0102     throw cms::Exception("SimTrackManager::saveTrackAndItsBranch") << " cannot handle hits for tracking";
0103   }
0104   trkH->setToBeSaved();
0105   int parent = trkH->parentID();
0106
0107   auto tk_itr =
0108       std::lower_bound(m_trackContainer.begin(), m_trackContainer.end(), parent, SimTrackManager::StrictWeakOrdering());
0109
0110   if (tk_itr != m_trackContainer.end() && (*tk_itr)->trackID() == parent) {
0111     saveTrackAndItsBranch(*tk_itr);
0112   }
0113 }
0114
0115 void SimTrackManager::storeTracks() {
0116   cleanTracksWithHistory();
0117
0118   // fill the map with the final mother-daughter relationship
0119   idsave.swap(ancestorList);
0120   std::stable_sort(idsave.begin(), idsave.end());
0121
0122   std::vector<std::pair<int, int> >().swap(ancestorList);
0123
0124   // to get a backward compatible order
0125   std::stable_sort(m_trackContainer.begin(), m_trackContainer.end(), trkIDLess());
0126
0127   // to reset the GenParticle ID of a SimTrack to its pre-LHCTransport value
0128   resetGenID();
0129
0130   reallyStoreTracks();
0131 }
0132
0133 void SimTrackManager::reallyStoreTracks() {
0134   // loop over the (now ordered) vector and really save the tracks
0135 #ifdef DebugLog
0136   edm::LogVerbatim("SimTrackManager") << "reallyStoreTracks() NtracksWithHistory= " << m_trackContainer.size();
0137 #endif
0138
0139   int nn = m_endPoints.size();
0140   for (auto& trkH : m_trackContainer) {
0141     // at this stage there is one vertex per track,
0142     // so the vertex id of track N is also N
0143     int iParentID = trkH->parentID();
0144     int ig = trkH->genParticleID();  // filled only for primary tracks
0145     bool isBackScatter = trkH->isFromBackScattering();
0146     bool isPrimary = trkH->isPrimary();
0147     int primaryGenPartId = trkH->getPrimaryID();  // filled if the G4Track had this info
0148     int ivertex = getOrCreateVertex(trkH, iParentID);
0149
0150     auto ptr = trkH;
0151     if (0 < iParentID) {
0152       for (auto& trk : m_trackContainer) {
0153         if (trk->trackID() == iParentID) {
0154           ptr = trk;
0155           break;
0156         }
0157       }
0158     }
0159     // Track at surface is the track at intersection point between tracker and calo
0160     // envelops if exist. If not exist the momentum is zero, position is the end of
0161     // the track
0162     const math::XYZVectorD& pm = ptr->momentum();
0163     math::XYZVectorD spos(0., 0., 0.);
0164     math::XYZTLorentzVectorD smom(0., 0., 0., 0.);
0165     int id = trkH->trackID();
0166     if (trkH->crossedBoundary()) {
0167       spos = trkH->getPositionAtBoundary();
0168       smom = trkH->getMomentumAtBoundary();
0169     } else {
0170       for (int i = 0; i < nn; ++i) {
0171         if (id == m_endPoints[i].first) {
0172           spos = m_endPoints[i].second;
0173           break;
0174         }
0175       }
0176     }
0177
0178     if (id >= static_cast<int>(PSimHit::k_tidOffset)) {
0179       throw cms::Exception("SimTrackManager::reallyStoreTracks")
0180           << " SimTrack ID " << id << " exceeds maximum allowed by PSimHit identifier" << PSimHit::k_tidOffset;
0181     }
0182     TmpSimTrack* g4simtrack =
0183         new TmpSimTrack(id, trkH->particleID(), trkH->momentum(), trkH->totalEnergy(), ivertex, ig, pm, spos, smom);
0184     g4simtrack->copyCrossedBoundaryVars(trkH);
0185     if (isBackScatter)
0186       g4simtrack->setFromBackScattering();
0187     if (isPrimary)
0188       g4simtrack->setIsPrimary();
0189     g4simtrack->setGenParticleID(primaryGenPartId);
0190     m_simEvent->addTrack(g4simtrack);
0191   }
0192 }
0193
0194 int SimTrackManager::getOrCreateVertex(TrackWithHistory* trkH, int iParentID) {
0195   int parent = -1;
0196   for (auto const& trk : m_trackContainer) {
0197     int id = trk->trackID();
0198     if (id == iParentID) {
0199       parent = id;
0200       break;
0201     }
0202   }
0203
0204   VertexMap::const_iterator iterator = m_vertexMap.find(parent);
0205   if (iterator != m_vertexMap.end()) {
0206     // loop over saved vertices
0207     for (auto const& xx : m_vertexMap[parent]) {
0208       if ((trkH->vertexPosition() - xx.second).Mag2() < r_limit2) {
0209         return xx.first;
0210       }
0211     }
0212   }
0213
0214   m_simEvent->addVertex(new TmpSimVertex(trkH->vertexPosition(), trkH->time(), parent, trkH->processType()));
0215   m_vertexMap[parent].push_back(VertexPosition(m_nVertices, trkH->vertexPosition()));
0216   ++m_nVertices;
0217   return (m_nVertices - 1);
0218 }
0219
0220 int SimTrackManager::idSavedTrack(int id) const {
0221   int idMother = id;
0222   if (id > 0) {
0223     unsigned int n = idsave.size();
0224     if (0 < n) {
0225       int jmax = n - 1;
0226       int j, id1;
0227
0228       // first loop forward
0229       bool notFound = true;
0230       for (j = 0; j <= jmax; ++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         return 0;
0244       }
0245
0246       // recursive loop
0247       do {
0248         notFound = true;
0249         // search ID scan backward
0250         for (j = jmax; j >= 0; --j) {
0251           if ((idsave[j]).first == idMother) {
0252             id1 = (idsave[j]).second;
0253             if (0 == id1 || id1 == idMother) {
0254               return id1;
0255             }
0256             jmax = j - 1;
0257             idMother = id1;
0258             notFound = false;
0259             break;
0260           }
0261         }
0262         if (notFound) {
0263           // ID not in the list of saved track - look into ancestors
0264           jmax = ancestorList.size() - 1;
0265           for (j = jmax; j >= 0; --j) {
0266             if ((ancestorList[j]).first == idMother) {
0267               idMother = (ancestorList[j]).second;
0268               return idMother;
0269             }
0270           }
0271           return 0;
0272         }
0273       } while (!notFound);
0274     }
0275   }
0276   return idMother;
0277 }
0278
0279 void SimTrackManager::fillMotherList() {
0280   if (!ancestorList.empty() && lastHist > ancestorList.size()) {
0281     lastHist = ancestorList.size();
0282     edm::LogError("SimTrackManager") << " SimTrackManager::fillMotherList track index corrupted";
0283   }
0284 #ifdef DebugLog
0285   edm::LogVerbatim("SimTrackManager") << "### SimTrackManager::fillMotherList: " << idsave.size()
0286                                       << " saved; ancestor: " << lastHist << "  " << ancestorList.size();
0287   for (unsigned int i = 0; i < idsave.size(); ++i) {
0288     edm::LogVerbatim("SimTrackManager") << " ISV: Track ID = " << (idsave[i]).first
0289                                         << " Mother ID = " << (idsave[i]).second;
0290   }
0291 #endif
0292   for (unsigned int n = lastHist; n < ancestorList.size(); ++n) {
0293     int theMotherId = idSavedTrack((ancestorList[n]).first);
0294     ancestorList[n].second = theMotherId;
0295 #ifdef DebugLog
0296     LogDebug("SimTrackManager") << "Track ID = " << (ancestorList[n]).first
0297                                 << " Mother ID = " << (ancestorList[n]).second;
0298 #endif
0299   }
0300
0301   lastHist = ancestorList.size();
0302   idsave.clear();
0303 }
0304
0305 void SimTrackManager::cleanTracksWithHistory() {
0306   if (m_trackContainer.empty() && idsave.empty()) {
0307     return;
0308   }
0309
0310 #ifdef DebugLog
0311   LogDebug("SimTrackManager") << "SimTrackManager::cleanTracksWithHistory has " << idsave.size()
0312                               << " mother-daughter relationships stored with lastTrack = " << lastTrack;
0313 #endif
0314
0315   if (lastTrack > 0 && lastTrack >= m_trackContainer.size()) {
0316     lastTrack = 0;
0317     edm::LogError("SimTrackManager") << " SimTrackManager::cleanTracksWithHistory track index corrupted";
0318   }
0319
0320   std::stable_sort(m_trackContainer.begin() + lastTrack, m_trackContainer.end(), trkIDLess());
0321   std::stable_sort(idsave.begin(), idsave.end());
0322
0323 #ifdef DebugLog
0324   LogDebug("SimTrackManager") << " SimTrackManager::cleanTracksWithHistory knows " << m_trackContainer.size()
0325                               << " tracks with history before branching";
0326   for (unsigned int it = 0; it < m_trackContainer.size(); it++) {
0327     LogDebug("SimTrackManager") << " 1 - Track in position " << it << " G4 track number "
0328                                 << m_trackContainer[it]->trackID() << " mother " << m_trackContainer[it]->parentID()
0329                                 << " status " << m_trackContainer[it]->saved();
0330   }
0331 #endif
0332
0333   for (auto const& t : m_trackContainer) {
0334     if (t->saved()) {
0335       saveTrackAndItsBranch(t);
0336     }
0337   }
0338   unsigned int num = lastTrack;
0339   for (unsigned int it = lastTrack; it < m_trackContainer.size(); ++it) {
0340     auto t = m_trackContainer[it];
0341     int g4ID = m_trackContainer[it]->trackID();
0342     if (t->saved()) {
0343       if (it > num) {
0344         m_trackContainer[num] = t;
0345       }
0346       ++num;
0347       for (auto& xx : idsave) {
0348         if (xx.first == g4ID) {
0349           xx.second = g4ID;
0350           break;
0351         }
0352       }
0353     } else {
0354       delete t;
0355     }
0356   }
0357
0358   m_trackContainer.resize(num);
0359
0360 #ifdef DebugLog
0361   LogDebug("SimTrackManager") << " AFTER CLEANING, I GET " << m_trackContainer.size()
0362                               << " tracks to be saved persistently";
0363   for (unsigned int it = 0; it < m_trackContainer.size(); ++it) {
0364     LogDebug("SimTrackManager") << " Track in position " << it << " G4 track number " << m_trackContainer[it]->trackID()
0365                                 << " mother " << m_trackContainer[it]->parentID() << " Status "
0366                                 << m_trackContainer[it]->saved() << " id " << m_trackContainer[it]->particleID()
0367                                 << " E(MeV)= " << m_trackContainer[it]->totalEnergy();
0368   }
0369 #endif
0370
0371   fillMotherList();
0372   lastTrack = m_trackContainer.size();
0373 }
0374
0375 void SimTrackManager::resetGenID() {
0376   if (theLHCTlink == nullptr)
0377     return;
0378
0379   for (auto const& trkH : m_trackContainer) {
0380     int genParticleID = trkH->genParticleID();
0381     if (genParticleID == -1) {
0382       continue;
0383     } else {
0384       for (auto const& xx : *theLHCTlink) {
0385         if (xx.afterHector() == genParticleID) {
0386           trkH->setGenParticleID(xx.beforeHector());
0387           continue;
0388         }
0389       }
0390     }
0391   }
0392
0393   theLHCTlink = nullptr;
0394 }
0395
0396 void SimTrackManager::ReportException(unsigned int id) const {
0397   throw cms::Exception("Unknown", "SimTrackManager::getTrackByID")
0398       << "Fail to get track " << id << " from SimTrackManager, container size= " << m_trackContainer.size();
0399 }