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     m_trackContainer.push_back(iTrack);
0077     const auto& v = track->GetStep()->GetPostStepPoint()->GetPosition();
0078     std::pair<int, math::XYZVectorD> p(iTrack->trackID(),
0079                                        math::XYZVectorD(v.x() * invcm, v.y() * invcm, v.z() * invcm));
0080     m_endPoints.push_back(p);
0081   }
0082   if (withAncestor) {
0083     std::pair<int, int> thisPair(iTrack->trackID(), 0);
0084     ancestorList.push_back(thisPair);
0085   }
0086 }
0087
0088 /// this saves a track and all its parents looping over the non ordered vector
0089 void SimTrackManager::saveTrackAndItsBranch(TrackWithHistory* trkWHist) {
0090   TrackWithHistory* trkH = trkWHist;
0091   if (trkH == nullptr) {
0092     edm::LogError("SimTrackManager") << " SimTrackManager::saveTrackAndItsBranch got 0 pointer ";
0093     throw cms::Exception("SimTrackManager::saveTrackAndItsBranch") << " cannot handle hits for tracking";
0094   }
0095   trkH->setToBeSaved();
0096   int parent = trkH->parentID();
0097
0098   auto tk_itr =
0099       std::lower_bound(m_trackContainer.begin(), m_trackContainer.end(), parent, SimTrackManager::StrictWeakOrdering());
0100
0101   if (tk_itr != m_trackContainer.end() && (*tk_itr)->trackID() == parent) {
0102     saveTrackAndItsBranch(*tk_itr);
0103   }
0104 }
0105
0106 void SimTrackManager::storeTracks() {
0107   cleanTracksWithHistory();
0108
0109   // fill the map with the final mother-daughter relationship
0110   idsave.swap(ancestorList);
0111   std::stable_sort(idsave.begin(), idsave.end());
0112
0113   std::vector<std::pair<int, int> >().swap(ancestorList);
0114
0115   // to get a backward compatible order
0116   std::stable_sort(m_trackContainer.begin(), m_trackContainer.end(), trkIDLess());
0117
0118   // to reset the GenParticle ID of a SimTrack to its pre-LHCTransport value
0119   resetGenID();
0120
0121   reallyStoreTracks();
0122 }
0123
0124 void SimTrackManager::reallyStoreTracks() {
0125   // loop over the (now ordered) vector and really save the tracks
0126 #ifdef DebugLog
0127   edm::LogVerbatim("SimTrackManager") << "reallyStoreTracks() NtracksWithHistory= " << m_trackContainer->size();
0128 #endif
0129
0130   int nn = m_endPoints.size();
0131   for (auto& trkH : m_trackContainer) {
0132     // at this stage there is one vertex per track,
0133     // so the vertex id of track N is also N
0134     int iParentID = trkH->parentID();
0135     int ig = trkH->genParticleID();
0136     int ivertex = getOrCreateVertex(trkH, iParentID);
0137
0138     auto ptr = trkH;
0139     if (0 < iParentID) {
0140       for (auto& trk : m_trackContainer) {
0141         if (trk->trackID() == iParentID) {
0142           ptr = trk;
0143           break;
0144         }
0145       }
0146     }
0147     // Track at surface is the track at intersection point between tracker and calo
0148     // envelops if exist. If not exist the momentum is zero, position is the end of
0149     // the track
0150     const math::XYZVectorD& pm = ptr->momentum();
0151     math::XYZVectorD spos(0., 0., 0.);
0152     math::XYZTLorentzVectorD smom(0., 0., 0., 0.);
0153     int id = trkH->trackID();
0154     if (trkH->crossedBoundary()) {
0155       spos = trkH->getPositionAtBoundary();
0156       smom = trkH->getMomentumAtBoundary();
0157     } else {
0158       for (int i = 0; i < nn; ++i) {
0159         if (id == m_endPoints[i].first) {
0160           spos = m_endPoints[i].second;
0161           break;
0162         }
0163       }
0164     }
0165
0166     if (id >= static_cast<int>(PSimHit::k_tidOffset)) {
0167       throw cms::Exception("SimTrackManager::reallyStoreTracks")
0168           << " SimTrack ID " << id << " exceeds maximum allowed by PSimHit identifier" << PSimHit::k_tidOffset;
0169     }
0170     TmpSimTrack* g4simtrack =
0171         new TmpSimTrack(id, trkH->particleID(), trkH->momentum(), trkH->totalEnergy(), ivertex, ig, pm, spos, smom);
0172     g4simtrack->copyCrossedBoundaryVars(trkH);
0173     m_simEvent->addTrack(g4simtrack);
0174   }
0175 }
0176
0177 int SimTrackManager::getOrCreateVertex(TrackWithHistory* trkH, int iParentID) {
0178   int parent = -1;
0179   for (auto const& trk : m_trackContainer) {
0180     int id = trk->trackID();
0181     if (id == iParentID) {
0182       parent = id;
0183       break;
0184     }
0185   }
0186
0187   VertexMap::const_iterator iterator = m_vertexMap.find(parent);
0188   if (iterator != m_vertexMap.end()) {
0189     // loop over saved vertices
0190     for (auto const& xx : m_vertexMap[parent]) {
0191       if ((trkH->vertexPosition() - xx.second).Mag2() < r_limit2) {
0192         return xx.first;
0193       }
0194     }
0195   }
0196
0197   m_simEvent->addVertex(new TmpSimVertex(trkH->vertexPosition(), trkH->time(), parent, trkH->processType()));
0198   m_vertexMap[parent].push_back(VertexPosition(m_nVertices, trkH->vertexPosition()));
0199   ++m_nVertices;
0200   return (m_nVertices - 1);
0201 }
0202
0203 int SimTrackManager::idSavedTrack(int id) const {
0204   int idMother = id;
0205   if (id > 0) {
0206     unsigned int n = idsave.size();
0207     if (0 < n) {
0208       int jmax = n - 1;
0209       int j, id1;
0210
0211       // first loop forward
0212       bool notFound = true;
0213       for (j = 0; j <= jmax; ++j) {
0214         if ((idsave[j]).first == idMother) {
0215           id1 = (idsave[j]).second;
0216           if (0 == id1 || id1 == idMother) {
0217             return id1;
0218           }
0219           jmax = j - 1;
0220           idMother = id1;
0221           notFound = false;
0222           break;
0223         }
0224       }
0225       if (notFound) {
0226         return 0;
0227       }
0228
0229       // recursive loop
0230       do {
0231         notFound = true;
0232         // search ID scan backward
0233         for (j = jmax; j >= 0; --j) {
0234           if ((idsave[j]).first == idMother) {
0235             id1 = (idsave[j]).second;
0236             if (0 == id1 || id1 == idMother) {
0237               return id1;
0238             }
0239             jmax = j - 1;
0240             idMother = id1;
0241             notFound = false;
0242             break;
0243           }
0244         }
0245         if (notFound) {
0246           // ID not in the list of saved track - look into ancestors
0247           jmax = ancestorList.size() - 1;
0248           for (j = jmax; j >= 0; --j) {
0249             if ((ancestorList[j]).first == idMother) {
0250               idMother = (ancestorList[j]).second;
0251               return idMother;
0252             }
0253           }
0254           return 0;
0255         }
0256       } while (!notFound);
0257     }
0258   }
0259   return idMother;
0260 }
0261
0262 void SimTrackManager::fillMotherList() {
0263   if (!ancestorList.empty() && lastHist > ancestorList.size()) {
0264     lastHist = ancestorList.size();
0265     edm::LogError("SimTrackManager") << " SimTrackManager::fillMotherList track index corrupted";
0266   }
0267 #ifdef DebugLog
0268   edm::LogVerbatim("SimTrackManager") << "### SimTrackManager::fillMotherList: " << idsave.size()
0269                                       << " saved; ancestor: " << lastHist << "  " << ancestorList.size();
0270   for (unsigned int i = 0; i < idsave.size(); ++i) {
0271     edm::LogVerbatim("SimTrackManager") << " ISV: Track ID = " << (idsave[i]).first
0272                                         << " Mother ID = " << (idsave[i]).second;
0273   }
0274 #endif
0275   for (unsigned int n = lastHist; n < ancestorList.size(); ++n) {
0276     int theMotherId = idSavedTrack((ancestorList[n]).first);
0277     ancestorList[n].second = theMotherId;
0278 #ifdef DebugLog
0279     LogDebug("SimTrackManager") << "Track ID = " << (ancestorList[n]).first
0280                                 << " Mother ID = " << (ancestorList[n]).second;
0281 #endif
0282   }
0283
0284   lastHist = ancestorList.size();
0285   idsave.clear();
0286 }
0287
0288 void SimTrackManager::cleanTracksWithHistory() {
0289   if (m_trackContainer.empty() && idsave.empty()) {
0290     return;
0291   }
0292
0293 #ifdef DebugLog
0294   LogDebug("SimTrackManager") << "SimTrackManager::cleanTracksWithHistory has " << idsave.size()
0295                               << " mother-daughter relationships stored with lastTrack = " << lastTrack;
0296 #endif
0297
0298   if (lastTrack > 0 && lastTrack >= m_trackContainer.size()) {
0299     lastTrack = 0;
0300     edm::LogError("SimTrackManager") << " SimTrackManager::cleanTracksWithHistory track index corrupted";
0301   }
0302
0303   std::stable_sort(m_trackContainer.begin() + lastTrack, m_trackContainer.end(), trkIDLess());
0304   std::stable_sort(idsave.begin(), idsave.end());
0305
0306 #ifdef DebugLog
0307   LogDebug("SimTrackManager") << " SimTrackManager::cleanTracksWithHistory knows " << m_trksForThisEvent->size()
0308                               << " tracks with history before branching";
0309   for (unsigned int it = 0; it < m_trackContainer.size(); it++) {
0310     LogDebug("SimTrackManager") << " 1 - Track in position " << it << " G4 track number "
0311                                 << m_trackContainer[it]->trackID() << " mother " << m_trackContainer[it]->parentID()
0312                                 << " status " << m_trackContainer[it]->saved();
0313   }
0314 #endif
0315
0316   for (auto const& t : m_trackContainer) {
0317     if (t->saved()) {
0318       saveTrackAndItsBranch(t);
0319     }
0320   }
0321   unsigned int num = lastTrack;
0322   for (unsigned int it = lastTrack; it < m_trackContainer.size(); ++it) {
0323     auto t = m_trackContainer[it];
0324     int g4ID = m_trackContainer[it]->trackID();
0325     if (t->saved()) {
0326       if (it > num) {
0327         m_trackContainer[num] = t;
0328       }
0329       ++num;
0330       for (auto& xx : idsave) {
0331         if (xx.first == g4ID) {
0332           xx.second = g4ID;
0333           break;
0334         }
0335       }
0336     } else {
0337       delete t;
0338     }
0339   }
0340
0341   m_trackContainer.resize(num);
0342
0343 #ifdef DebugLog
0344   LogDebug("SimTrackManager") << " AFTER CLEANING, I GET " << m_trackContainer.size()
0345                               << " tracks to be saved persistently";
0346   for (unsigned int it = 0; it < m_trackContainer.size(); ++it) {
0347     LogDebug("SimTrackManager") << " Track in position " << it << " G4 track number " << m_trackContainer[it]->trackID()
0348                                 << " mother " << m_trackContainer[it]->parentID() << " Status "
0349                                 << m_trackContainer[it]->saved() << " id " << m_trackContainer[it]->particleID()
0350                                 << " E(MeV)= " << m_trackContainer[it]->totalEnergy();
0351   }
0352 #endif
0353
0354   fillMotherList();
0355   lastTrack = m_trackContainer.size();
0356 }
0357
0358 void SimTrackManager::resetGenID() {
0359   if (theLHCTlink == nullptr)
0360     return;
0361
0362   for (auto const& trkH : m_trackContainer) {
0363     int genParticleID = trkH->genParticleID();
0364     if (genParticleID == -1) {
0365       continue;
0366     } else {
0367       for (auto const& xx : *theLHCTlink) {
0368         if (xx.afterHector() == genParticleID) {
0369           trkH->setGenParticleID(xx.beforeHector());
0370           continue;
0371         }
0372       }
0373     }
0374   }
0375
0376   theLHCTlink = nullptr;
0377 }
0378
0379 void SimTrackManager::ReportException(unsigned int id) const {
0380   throw cms::Exception("Unknown", "SimTrackManager::getTrackByID")
0381       << "Fail to get track " << id << " from SimTrackManager, container size= " << m_trackContainer.size();
0382 }