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/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 //#define DebugLog
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 // member functions
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 /// this saves a track and all its parents looping over the non ordered vector
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   // fill the map with the final mother-daughter relationship
0093   idsave.swap(ancestorList);
0094   stable_sort(idsave.begin(), idsave.end());
0095
0096   std::vector<std::pair<int, int> >().swap(ancestorList);
0097
0098   // to get a backward compatible order
0099   stable_sort(m_trksForThisEvent->begin(), m_trksForThisEvent->end(), trkIDLess());
0100
0101   // to reset the GenParticle ID of a SimTrack to its pre-LHCTransport value
0102   resetGenID();
0103
0104   reallyStoreTracks(simEvent);
0105 }
0106
0107 void SimTrackManager::reallyStoreTracks(G4SimEvent* simEvent) {
0108   // loop over the (now ordered) vector and really save the tracks
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     // at this stage there is one vertex per track,
0116     // so the vertex id of track N is also N
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     // loop over saved vertices
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       // first loop forward
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       // recursive loop
0219       do {
0220         notFound = true;
0221         // search ID scan backward
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           // ID not in the list of saved track - look into ancestors
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   edm::LogVerbatim("SimTrackManager") << "### SimTrackManager::fillMotherList: "
0258         << idsave.size() << " saved; ancestor: " << lastHist
0259         << "  " << ancestorList.size() << std::endl;
0260   for (unsigned int i = 0; i< idsave.size(); ++i) {
0261     edm::LogVerbatim("SimTrackManager")  << " ISV: Track ID = " << (idsave[i]).first
0262            << " Mother ID = " << (idsave[i]).second << std::endl;
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     std::cout  << " ANC: Track ID = " << (ancestorList[n]).first
0270            << " Mother ID = " << (ancestorList[n]).second << std::endl;
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 }