diff --git a/SimDataFormats/Track/interface/SimTrack.h b/SimDataFormats/Track/interface/SimTrack.h index dc5dce6ca3646..19f73ea6b385a 100644 --- a/SimDataFormats/Track/interface/SimTrack.h +++ b/SimDataFormats/Track/interface/SimTrack.h @@ -34,8 +34,9 @@ class SimTrack : public CoreSimTrack { bool noVertex() const { return ivert == -1; } /// index of the corresponding Generator particle in the Event container (-1 if no Genpart) - int genpartIndex() const { return igenpart; } - bool noGenpart() const { return igenpart == -1; } + bool isPrimary() const { return (trackInfo_ >> 1) & 1; } + int genpartIndex() const { return isPrimary() ? igenpart : -1; } + bool noGenpart() const { return isPrimary() ? igenpart == -1 : true; } const math::XYZVectorD& trackerSurfacePosition() const { return tkposition; } @@ -51,16 +52,25 @@ class SimTrack : public CoreSimTrack { int idAtBoundary, math::XYZTLorentzVectorF positionAtBoundary, math::XYZTLorentzVectorF momentumAtBoundary) { - crossedBoundary_ = crossedBoundary; + if (crossedBoundary) + trackInfo_ |= (1 << 2); idAtBoundary_ = idAtBoundary; positionAtBoundary_ = positionAtBoundary; momentumAtBoundary_ = momentumAtBoundary; } - bool crossedBoundary() const { return crossedBoundary_; } + bool crossedBoundary() const { return (trackInfo_ >> 2) & 1; } const math::XYZTLorentzVectorF& getPositionAtBoundary() const { return positionAtBoundary_; } const math::XYZTLorentzVectorF& getMomentumAtBoundary() const { return momentumAtBoundary_; } int getIDAtBoundary() const { return idAtBoundary_; } + bool isFromBackScattering() const { return trackInfo_ & 1; } + void setFromBackScattering() { trackInfo_ |= 1; } + + void setIsPrimary() { trackInfo_ |= (1 << 1); } + void setGenParticleID(const int idx) { igenpart = idx; } + int getPrimaryID() const { return igenpart; } + uint8_t getTrackInfo() const { return trackInfo_; } + private: int ivert; int igenpart; @@ -68,10 +78,14 @@ class SimTrack : public CoreSimTrack { math::XYZVectorD tkposition; math::XYZTLorentzVectorD tkmomentum; - bool crossedBoundary_; int idAtBoundary_; math::XYZTLorentzVectorF positionAtBoundary_; math::XYZTLorentzVectorF momentumAtBoundary_; + uint8_t trackInfo_; + // explanation of trackInfo bits: + // 00000001 = simTrack is from backscattering + // 00000010 = simTrack is of a primary particle + // 00000100 = simTrack crossed the boundary }; #include diff --git a/SimDataFormats/Track/src/SimTrack.cc b/SimDataFormats/Track/src/SimTrack.cc index 52b3441e719b6..6126d4d6fc6b0 100644 --- a/SimDataFormats/Track/src/SimTrack.cc +++ b/SimDataFormats/Track/src/SimTrack.cc @@ -1,12 +1,12 @@ #include "SimDataFormats/Track/interface/SimTrack.h" -SimTrack::SimTrack() : ivert(-1), igenpart(-1), crossedBoundary_(false) {} +SimTrack::SimTrack() : ivert(-1), igenpart(-1), trackInfo_(0) {} SimTrack::SimTrack(int ipart, const math::XYZTLorentzVectorD& p) - : Core(ipart, p), ivert(-1), igenpart(-1), crossedBoundary_(false) {} + : Core(ipart, p), ivert(-1), igenpart(-1), trackInfo_(0) {} SimTrack::SimTrack(int ipart, const math::XYZTLorentzVectorD& p, int iv, int ig) - : Core(ipart, p), ivert(iv), igenpart(ig), crossedBoundary_(false) {} + : Core(ipart, p), ivert(iv), igenpart(ig), trackInfo_(0) {} SimTrack::SimTrack(int ipart, const math::XYZTLorentzVectorD& p, @@ -14,9 +14,9 @@ SimTrack::SimTrack(int ipart, int ig, const math::XYZVectorD& tkp, const math::XYZTLorentzVectorD& tkm) - : Core(ipart, p), ivert(iv), igenpart(ig), tkposition(tkp), tkmomentum(tkm), crossedBoundary_(false) {} + : Core(ipart, p), ivert(iv), igenpart(ig), tkposition(tkp), tkmomentum(tkm), trackInfo_(0) {} -SimTrack::SimTrack(const CoreSimTrack& t, int iv, int ig) : Core(t), ivert(iv), igenpart(ig), crossedBoundary_(false) {} +SimTrack::SimTrack(const CoreSimTrack& t, int iv, int ig) : Core(t), ivert(iv), igenpart(ig), trackInfo_(0) {} std::ostream& operator<<(std::ostream& o, const SimTrack& t) { return o << (SimTrack::Core)(t) << ", " << t.vertIndex() << ", " << t.genpartIndex(); diff --git a/SimDataFormats/Track/src/classes_def.xml b/SimDataFormats/Track/src/classes_def.xml index d76801ad762c7..09c6a03d5247e 100644 --- a/SimDataFormats/Track/src/classes_def.xml +++ b/SimDataFormats/Track/src/classes_def.xml @@ -4,10 +4,27 @@ - + + + + type(), newObj->momentum(), newObj->vertIndex(), newObj->genpartIndex(), newObj->trackerSurfacePosition(), newObj->trackerSurfaceMomentum()); + tmp.setTrackId(newObj->trackId()); + tmp.setEventId(newObj->eventId()); + tmp.setCrossedBoundaryVars( + newObj->crossedBoundary(), newObj->getIDAtBoundary(), newObj->getPositionAtBoundary(), newObj->getMomentumAtBoundary()); + if (newObj->isFromBackScattering()) { + tmp.setFromBackScattering(); + } + if (newObj->genpartIndex() != -1) { + tmp.setIsPrimary(); + } + *newObj=tmp; + ]]> + diff --git a/SimG4Core/Application/plugins/OscarMTProducer.cc b/SimG4Core/Application/plugins/OscarMTProducer.cc index 39e70c5823255..f465329c8ac45 100644 --- a/SimG4Core/Application/plugins/OscarMTProducer.cc +++ b/SimG4Core/Application/plugins/OscarMTProducer.cc @@ -251,7 +251,7 @@ void OscarMTProducer::produce(edm::Event& e, const edm::EventSetup& es) { if (0 < m_verbose) { edm::LogVerbatim("SimG4CoreApplication") - << "Produced " << p2->size() << " SimVertecies: position(cm), time(s), parentID, vertexID, processType"; + << "Produced " << p2->size() << " SimVertices: position(cm), time(s), parentID, vertexID, processType"; if (1 < m_verbose) { int nn = p2->size(); for (int i = 0; i < nn; ++i) { @@ -260,12 +260,15 @@ void OscarMTProducer::produce(edm::Event& e, const edm::EventSetup& es) { } edm::LogVerbatim("SimG4CoreApplication") << "Produced " << p1->size() - << " SimTracks: pdg, 4-momentum(GeV), vertexID, mcTruthID, flagBoundary, trackID at boundary"; + << " SimTracks: G4 Id, pdg, 4-momentum(GeV), vertexID, mcTruthID, crossedBoundary -> trackID at boundary, from " + "backscattering, isPrimary -> getPrimary"; if (1 < m_verbose) { int nn = p1->size(); for (int i = 0; i < nn; ++i) { - edm::LogVerbatim("Track") << " " << i << ". " << (*p1)[i] << " " << (*p1)[i].crossedBoundary() << " " - << (*p1)[i].getIDAtBoundary(); + edm::LogVerbatim("Track") << " " << i << ". " << (*p1)[i].trackId() << ", " << (*p1)[i] << ", " + << (*p1)[i].crossedBoundary() << "-> " << (*p1)[i].getIDAtBoundary() << ", " + << (*p1)[i].isFromBackScattering() << ", " << (*p1)[i].isPrimary() << "-> " + << (*p1)[i].getPrimaryID(); } } } diff --git a/SimG4Core/Notification/interface/TmpSimTrack.h b/SimG4Core/Notification/interface/TmpSimTrack.h index 4ecbe7a47a4f8..523a2825f195c 100644 --- a/SimG4Core/Notification/interface/TmpSimTrack.h +++ b/SimG4Core/Notification/interface/TmpSimTrack.h @@ -40,7 +40,7 @@ class TmpSimTrack { const math::XYZVectorD& momentum() const { return ip_; } double energy() const { return ie_; } int ivert() const { return ivert_; } - int igenpart() const { return igenpart_; } + int igenpart() const { return isPrimary_ ? igenpart_ : -1; } // parent momentum at interaction const math::XYZVectorD& parentMomentum() const { return parentMomentum_; } // Information at level of tracker surface @@ -62,6 +62,12 @@ class TmpSimTrack { const math::XYZTLorentzVectorF& getPositionAtBoundary() const { return positionAtBoundary_; } const math::XYZTLorentzVectorF& getMomentumAtBoundary() const { return momentumAtBoundary_; } int getIDAtBoundary() const { return idAtBoundary_; } + bool isFromBackScattering() const { return isFromBackScattering_; } + void setFromBackScattering() { isFromBackScattering_ = true; } + void setGenParticleID(int i) { igenpart_ = i; } + bool isPrimary() const { return isPrimary_; } + void setIsPrimary() { isPrimary_ = true; } + int getPrimaryID() const { return igenpart_; } private: int id_; @@ -74,7 +80,9 @@ class TmpSimTrack { math::XYZVectorD parentMomentum_{math::XYZVectorD(0., 0., 0.)}; math::XYZVectorD tkSurfacePosition_{math::XYZVectorD(0., 0., 0.)}; math::XYZTLorentzVectorD tkSurfaceMomentum_{math::XYZTLorentzVectorD(0., 0., 0., 0.)}; + bool isFromBackScattering_{false}; bool crossedBoundary_{false}; + bool isPrimary_{false}; int idAtBoundary_{-1}; math::XYZTLorentzVectorF positionAtBoundary_{math::XYZTLorentzVectorF(0.f, 0.f, 0.f, 0.f)}; math::XYZTLorentzVectorF momentumAtBoundary_{math::XYZTLorentzVectorF(0.f, 0.f, 0.f, 0.f)}; diff --git a/SimG4Core/Notification/interface/TrackWithHistory.h b/SimG4Core/Notification/interface/TrackWithHistory.h index 4c272cf301c5f..f82bce37b105c 100644 --- a/SimG4Core/Notification/interface/TrackWithHistory.h +++ b/SimG4Core/Notification/interface/TrackWithHistory.h @@ -14,7 +14,7 @@ class G4PrimaryParticle; class TrackWithHistory { public: - /** The constructor is called at time, + /** The constructor is called at time, * when some of the information may not available yet. */ TrackWithHistory(const G4Track *g4track, int pID); @@ -27,7 +27,7 @@ class TrackWithHistory { int trackID() const { return trackID_; } int particleID() const { return pdgID_; } int parentID() const { return parentID_; } - int genParticleID() const { return genParticleID_; } + int genParticleID() const { return isPrimary_ ? genParticleID_ : -1; } int vertexID() const { return vertexID_; } int processType() const { return procType_; } int getIDAtBoundary() const { return idAtBoundary_; } @@ -67,6 +67,11 @@ class TrackWithHistory { tkSurfacePosition_ = pos; tkSurfaceMomentum_ = mom; } + bool isFromBackScattering() const { return isFromBackScattering_; } + void setFromBackScattering() { isFromBackScattering_ = true; } + bool isPrimary() const { return isPrimary_; } + void setIsPrimary() { isPrimary_ = true; } + int getPrimaryID() const { return genParticleID_; } private: int trackID_; @@ -88,6 +93,8 @@ class TrackWithHistory { bool storeTrack_{false}; bool saved_{false}; bool crossedBoundary_{false}; + bool isFromBackScattering_{false}; + bool isPrimary_{false}; }; extern G4ThreadLocal G4Allocator *fpTrackWithHistoryAllocator; diff --git a/SimG4Core/Notification/src/SimTrackManager.cc b/SimG4Core/Notification/src/SimTrackManager.cc index d8e05af69df9f..db4ef4006ca23 100644 --- a/SimG4Core/Notification/src/SimTrackManager.cc +++ b/SimG4Core/Notification/src/SimTrackManager.cc @@ -73,6 +73,15 @@ void SimTrackManager::addTrack(TrackWithHistory* iTrack, const G4Track* track, b std::pair thePair(iTrack->trackID(), iTrack->parentID()); idsave.push_back(thePair); if (inHistory) { + auto info = static_cast(track->GetUserInformation()); + if (info->isInTrkFromBackscattering()) + iTrack->setFromBackScattering(); + // set there for the *non-primary* tracks the genParticle ID associated with the G4Track + // for the primaries this is done in the TrackWithHistory constructor. + // In the constructor of TrackWithHistory the info isPrimary is saved and used + // to give -1 if the track is not a primary. + if (not iTrack->isPrimary()) + iTrack->setGenParticleID(info->mcTruthID()); m_trackContainer.push_back(iTrack); const auto& v = track->GetStep()->GetPostStepPoint()->GetPosition(); std::pair p(iTrack->trackID(), @@ -124,7 +133,7 @@ void SimTrackManager::storeTracks() { void SimTrackManager::reallyStoreTracks() { // loop over the (now ordered) vector and really save the tracks #ifdef DebugLog - edm::LogVerbatim("SimTrackManager") << "reallyStoreTracks() NtracksWithHistory= " << m_trackContainer->size(); + edm::LogVerbatim("SimTrackManager") << "reallyStoreTracks() NtracksWithHistory= " << m_trackContainer.size(); #endif int nn = m_endPoints.size(); @@ -132,7 +141,10 @@ void SimTrackManager::reallyStoreTracks() { // at this stage there is one vertex per track, // so the vertex id of track N is also N int iParentID = trkH->parentID(); - int ig = trkH->genParticleID(); + int ig = trkH->genParticleID(); // filled only for primary tracks + bool isBackScatter = trkH->isFromBackScattering(); + bool isPrimary = trkH->isPrimary(); + int primaryGenPartId = trkH->getPrimaryID(); // filled if the G4Track had this info int ivertex = getOrCreateVertex(trkH, iParentID); auto ptr = trkH; @@ -170,6 +182,11 @@ void SimTrackManager::reallyStoreTracks() { TmpSimTrack* g4simtrack = new TmpSimTrack(id, trkH->particleID(), trkH->momentum(), trkH->totalEnergy(), ivertex, ig, pm, spos, smom); g4simtrack->copyCrossedBoundaryVars(trkH); + if (isBackScatter) + g4simtrack->setFromBackScattering(); + if (isPrimary) + g4simtrack->setIsPrimary(); + g4simtrack->setGenParticleID(primaryGenPartId); m_simEvent->addTrack(g4simtrack); } } @@ -304,7 +321,7 @@ void SimTrackManager::cleanTracksWithHistory() { std::stable_sort(idsave.begin(), idsave.end()); #ifdef DebugLog - LogDebug("SimTrackManager") << " SimTrackManager::cleanTracksWithHistory knows " << m_trksForThisEvent->size() + LogDebug("SimTrackManager") << " SimTrackManager::cleanTracksWithHistory knows " << m_trackContainer.size() << " tracks with history before branching"; for (unsigned int it = 0; it < m_trackContainer.size(); it++) { LogDebug("SimTrackManager") << " 1 - Track in position " << it << " G4 track number " diff --git a/SimG4Core/Notification/src/TmpSimEvent.cc b/SimG4Core/Notification/src/TmpSimEvent.cc index 5d82d4a228adf..0d476825291ff 100644 --- a/SimG4Core/Notification/src/TmpSimEvent.cc +++ b/SimG4Core/Notification/src/TmpSimEvent.cc @@ -35,15 +35,26 @@ void TmpSimEvent::load(edm::SimTrackContainer& c) const { int iv = trk->ivert(); int ig = trk->igenpart(); int id = trk->id(); + bool isBackScatter = trk->isFromBackScattering(); + bool isPrimary = trk->isPrimary(); + int primaryGenPartId = trk->getPrimaryID(); // filled if the G4Track had this info // ip = particle ID as PDG - // pp = 4-momentum in GeV + // p = 4-momentum in GeV // iv = corresponding TmpSimVertex index - // ig = corresponding GenParticle index + // ig = corresponding GenParticle index only if is a primary, otherwise is -1 + // primaryGenPartId = corresponding GenParticle index also if not primary + // id = corresponding g4Track Id SimTrack t = SimTrack(ip, p, iv, ig, trk->trackerSurfacePosition(), trk->trackerSurfaceMomentum()); t.setTrackId(id); t.setEventId(EncodedEventId(0)); t.setCrossedBoundaryVars( trk->crossedBoundary(), trk->getIDAtBoundary(), trk->getPositionAtBoundary(), trk->getMomentumAtBoundary()); + if (isBackScatter) + t.setFromBackScattering(); + if (isPrimary) + t.setIsPrimary(); + else + t.setGenParticleID(primaryGenPartId); c.push_back(t); } std::stable_sort(c.begin(), c.end(), IdSort()); diff --git a/SimG4Core/Notification/src/TrackWithHistory.cc b/SimG4Core/Notification/src/TrackWithHistory.cc index ea6a843d19d51..060511bb154d1 100644 --- a/SimG4Core/Notification/src/TrackWithHistory.cc +++ b/SimG4Core/Notification/src/TrackWithHistory.cc @@ -29,10 +29,15 @@ TrackWithHistory::TrackWithHistory(const G4Track* g4trk, int pID) { TrackInformation* trkinfo = static_cast(g4trk->GetUserInformation()); storeTrack_ = trkinfo->storeTrack(); auto vgprimary = g4trk->GetDynamicParticle()->GetPrimaryParticle(); + // GetPrimaryParticle() returns the pointer to the corresponding G4PrimaryParticle object + // if this particle is a primary particle OR is defined as a + // pre-assigned decay product. Otherwise return nullptr. + // Therefore here the genParticleID_ is -1 for not primary particles if (vgprimary != nullptr) { auto priminfo = static_cast(vgprimary->GetUserInformation()); if (nullptr != priminfo) { genParticleID_ = priminfo->id(); + isPrimary_ = true; } } // V.I. weight is computed in the same way as before @@ -53,6 +58,7 @@ TrackWithHistory::TrackWithHistory(const G4PrimaryParticle* ptr, int trackID, co auto priminfo = static_cast(ptr->GetUserInformation()); if (nullptr != priminfo) { genParticleID_ = priminfo->id(); + isPrimary_ = true; } weight_ = 10000. * genParticleID_; }