diff --git a/RecoMTD/TimingIDTools/plugins/MVATrainingNtuple.cc b/RecoMTD/TimingIDTools/plugins/MVATrainingNtuple.cc index ccfbde216417b..aee7654e3fadc 100644 --- a/RecoMTD/TimingIDTools/plugins/MVATrainingNtuple.cc +++ b/RecoMTD/TimingIDTools/plugins/MVATrainingNtuple.cc @@ -709,9 +709,9 @@ void MVATrainingNtuple::analyze(const edm::Event& iEvent, const edm::EventSetup& }); // Check if TP has direct or other sim cluster for BTL for (const auto& simClusterRef : simClustersRefs) { - if (simClusterRef->trackIdOffset() == 0) { + if (simClusterRef->hitProdType() == 0) { isTPmtdDirectETL = true; - } else if (simClusterRef->trackIdOffset() != 0) { + } else if (simClusterRef->hitProdType() != 0) { isTPmtdOtherETL = true; } } @@ -724,9 +724,9 @@ void MVATrainingNtuple::analyze(const edm::Event& iEvent, const edm::EventSetup& }); // Check if TP has direct or other sim cluster for BTL for (const auto& simClusterRef : simClustersRefs) { - if (simClusterRef->trackIdOffset() == 0) { + if (simClusterRef->hitProdType() == 0) { isTPmtdDirectBTL = true; - } else if (simClusterRef->trackIdOffset() != 0) { + } else if (simClusterRef->hitProdType() != 0) { isTPmtdOtherBTL = true; } } diff --git a/SimDataFormats/CaloAnalysis/interface/MtdSimCluster.h b/SimDataFormats/CaloAnalysis/interface/MtdSimCluster.h index b7db1f4afaf60..a7502a65d6d90 100644 --- a/SimDataFormats/CaloAnalysis/interface/MtdSimCluster.h +++ b/SimDataFormats/CaloAnalysis/interface/MtdSimCluster.h @@ -116,9 +116,9 @@ class MtdSimCluster : public SimCluster { ++nsimhits_; } - void setTrackIdOffset(unsigned int offset) { idOffset_ = offset; } + void setHitProdType(unsigned int offset) { idOffset_ = offset; } - unsigned int trackIdOffset() const { return idOffset_; } + unsigned int hitProdType() const { return idOffset_; } protected: std::vector mtdHits_; diff --git a/SimDataFormats/CaloAnalysis/src/MtdSimLayerCluster.cc b/SimDataFormats/CaloAnalysis/src/MtdSimLayerCluster.cc index 201136111fca8..957711f1e822d 100644 --- a/SimDataFormats/CaloAnalysis/src/MtdSimLayerCluster.cc +++ b/SimDataFormats/CaloAnalysis/src/MtdSimLayerCluster.cc @@ -27,7 +27,7 @@ MtdSimLayerCluster::~MtdSimLayerCluster() {} std::ostream &operator<<(std::ostream &s, MtdSimLayerCluster const &tp) { s << "CP momentum, q, ID, & Event #: " << tp.p4() << " " << tp.charge() << " " << tp.pdgId() << " " << tp.eventId().bunchCrossing() << "." << tp.eventId().event() << std::endl; - s << " Offset " << tp.trackIdOffset() << " " + s << " Offset " << tp.hitProdType() << " " << " LC time " << tp.simLCTime() << " LC energy " << tp.simLCEnergy() << std::endl; for (MtdSimLayerCluster::genp_iterator hepT = tp.genParticle_begin(); hepT != tp.genParticle_end(); ++hepT) { diff --git a/SimDataFormats/TrackingHit/interface/PSimHit.h b/SimDataFormats/TrackingHit/interface/PSimHit.h index 375c7e0419cd3..d92aa15ad5994 100644 --- a/SimDataFormats/TrackingHit/interface/PSimHit.h +++ b/SimDataFormats/TrackingHit/interface/PSimHit.h @@ -14,8 +14,6 @@ class TrackingSlaveSD; // for friend declaration only class PSimHit { public: - static constexpr unsigned int k_tidOffset = 200000000; - PSimHit() : theDetUnitId(0) {} PSimHit(const Local3DPoint& entry, @@ -107,15 +105,6 @@ class PSimHit { */ unsigned int trackId() const { return theTrackId; } - /** In case te SimTrack ID is incremented by the k_tidOffset for hit category definition, this - * methods returns the original theTrackId value directly. - */ - unsigned int originalTrackId() const { return theTrackId % k_tidOffset; } - - unsigned int offsetTrackId() const { return theTrackId / k_tidOffset; } - - static unsigned int addTrackIdOffset(unsigned int tId, unsigned int offset) { return offset * k_tidOffset + tId; } - EncodedEventId eventId() const { return theEventId; } void setEventId(EncodedEventId e) { theEventId = e; } @@ -128,10 +117,20 @@ class PSimHit { * value with special significance is zero (for "undefined"), so zero should * not be the ID of any process. */ - unsigned short processType() const { return theProcessType; } + + // use 9 bits (up to 511) for process id, reserve the rest for hit production mechanism id + // 7 bits field available in PSimHit processType, i.e. up 127, to store processes + unsigned short processType() const { return theProcessType & kProcidMask; } + + unsigned short hitProdType() const { return (theProcessType >> kHitidShift) & kHitidMask; } + void setHitProdType(unsigned int hitId) { theProcessType |= hitId << kHitidShift; } void setTof(float tof) { theTof = tof; } + static constexpr unsigned int kProcidMask = 0x1FF; + static constexpr unsigned int kHitidMask = 0x7F; + static constexpr unsigned int kHitidShift = 9; + protected: // properties Local3DPoint theEntryPoint; // position at entry diff --git a/SimDataFormats/TrackingHit/interface/SimHitCategory.h b/SimDataFormats/TrackingHit/interface/SimHitCategory.h new file mode 100644 index 0000000000000..d62e3f34401a7 --- /dev/null +++ b/SimDataFormats/TrackingHit/interface/SimHitCategory.h @@ -0,0 +1,22 @@ +#ifndef SimHitCategory_h +#define SimHitCategory_h + +namespace SimHitCategory { + + // Identification code of sim hit production mechanism, subdetector dependent + + // MTD + + static constexpr unsigned int nCategoriesMTD = 5; + + static constexpr unsigned int prodTypeMTD[nCategoriesMTD] = { + 0, // direct hit from particle coming from tracker + 1, // BTL hit from secondary particle not saved in history + 2, // BTL hit from identified looper + 3, // BTL hit from back-scattering from CALO volume + 4 // ETL hit entering from a rear face of disks + }; + +}; // namespace SimHitCategory + +#endif diff --git a/SimDataFormats/TrackingHit/src/classes_def.xml b/SimDataFormats/TrackingHit/src/classes_def.xml index 9c815102df851..dcae308894bd5 100644 --- a/SimDataFormats/TrackingHit/src/classes_def.xml +++ b/SimDataFormats/TrackingHit/src/classes_def.xml @@ -1,6 +1,17 @@ - + + + + > 25) == 0x31) {theTrackId = (onfile.theTrackId % 200000000);} + ]]> + + + > 25) == 0x31) {unsigned short tmp = onfile.theProcessType; tmp |= (onfile.theTrackId / 200000000) << 9 ; theProcessType = tmp;} + ]]> + diff --git a/SimDataFormats/TrackingHit/test/BuildFile.xml b/SimDataFormats/TrackingHit/test/BuildFile.xml new file mode 100644 index 0000000000000..776344822a906 --- /dev/null +++ b/SimDataFormats/TrackingHit/test/BuildFile.xml @@ -0,0 +1,4 @@ + + + + diff --git a/SimDataFormats/TrackingHit/test/testPSimHit.cc b/SimDataFormats/TrackingHit/test/testPSimHit.cc new file mode 100644 index 0000000000000..986837c0cb149 --- /dev/null +++ b/SimDataFormats/TrackingHit/test/testPSimHit.cc @@ -0,0 +1,24 @@ +#include "SimDataFormats/TrackingHit/interface/PSimHit.h" + +#include +#include +#include + +int main() { + Local3DPoint dummy(0, 0, 0); + + for (int procId = 1; procId < 404; procId++) { + for (unsigned short itype = 0; itype < 8; itype++) { + PSimHit testHit(dummy, dummy, 0., 0., 0., 11, 0, 0, 0., 0., procId); + testHit.setHitProdType(itype); + std::cout << " hit procId = " << std::fixed << std::setw(8) << procId << " sim procId = " << std::setw(8) + << testHit.processType() << " hit type = " << std::setw(8) << itype << " sim type = " << std::setw(8) + << testHit.hitProdType() << std::endl; + + assert(procId == testHit.processType()); + assert(itype == testHit.hitProdType()); + } + } + + return 0; +} diff --git a/SimG4CMS/Forward/interface/BscG4Hit.h b/SimG4CMS/Forward/interface/BscG4Hit.h index 4aae1b75e3310..2e76af57b4867 100644 --- a/SimG4CMS/Forward/interface/BscG4Hit.h +++ b/SimG4CMS/Forward/interface/BscG4Hit.h @@ -14,6 +14,8 @@ #include #include +#include "SimDataFormats/TrackingHit/interface/PSimHit.h" + class BscG4Hit : public G4VHit { public: BscG4Hit(); @@ -92,6 +94,8 @@ class BscG4Hit : public G4VHit { int getParentId() const { return theParentId; }; int getProcessId() const { return theProcessId; }; + int getProdType() const { return (theProcessId >> PSimHit::kHitidShift) & PSimHit::kHitidMask; } + void setHitProdType(unsigned int hitId) { theProcessId |= hitId << PSimHit::kHitidShift; } float getVx() const { return theVx; }; float getVy() const { return theVy; }; float getVz() const { return theVz; }; diff --git a/SimG4CMS/Forward/interface/MtdHitCategory.h b/SimG4CMS/Forward/interface/MtdHitCategory.h deleted file mode 100644 index 29c4f02601249..0000000000000 --- a/SimG4CMS/Forward/interface/MtdHitCategory.h +++ /dev/null @@ -1,15 +0,0 @@ -#ifndef SimG4CMSForward_MtdHitCategory_h -#define SimG4CMSForward_MtdHitCategory_h - -#include - -namespace MtdHitCategory { - static constexpr unsigned int k_idsecOffset = 1; - static constexpr unsigned int k_idloopOffset = 2; - static constexpr unsigned int k_idFromCaloOffset = 3; - static constexpr unsigned int k_idETLfromBack = 4; - static constexpr unsigned int n_categories = - std::max({k_idsecOffset, k_idloopOffset, k_idFromCaloOffset, k_idETLfromBack}); -}; // namespace MtdHitCategory - -#endif diff --git a/SimG4CMS/Forward/interface/MtdSD.h b/SimG4CMS/Forward/interface/MtdSD.h index 19d651fb84996..b82f9d3a16afa 100644 --- a/SimG4CMS/Forward/interface/MtdSD.h +++ b/SimG4CMS/Forward/interface/MtdSD.h @@ -2,7 +2,7 @@ #define SimG4CMSForward_MtdSD_h #include "SimG4CMS/Forward/interface/TimingSD.h" -#include "SimG4CMS/Forward/interface/MtdHitCategory.h" +#include "SimDataFormats/TrackingHit/interface/SimHitCategory.h" #include "FWCore/Framework/interface/Frameworkfwd.h" #include "FWCore/ParameterSet/interface/ParameterSetfwd.h" diff --git a/SimG4CMS/Forward/src/BSCG4Hit.cc b/SimG4CMS/Forward/src/BSCG4Hit.cc index c8b0f74e162c4..f6f6163d9b955 100644 --- a/SimG4CMS/Forward/src/BSCG4Hit.cc +++ b/SimG4CMS/Forward/src/BSCG4Hit.cc @@ -5,6 +5,7 @@ /////////////////////////////////////////////////////////////////////////////// #include "FWCore/MessageLogger/interface/MessageLogger.h" #include "SimG4CMS/Forward/interface/BscG4Hit.h" +#include "SimDataFormats/TrackingHit/interface/SimHitCategory.h" #include BscG4Hit::BscG4Hit() : entry(0., 0., 0.), entrylp(0., 0., 0.), exitlp(0., 0., 0.) { @@ -141,8 +142,8 @@ std::ostream& operator<<(std::ostream& os, const BscG4Hit& hit) { << " EnergyLoss = " << hit.getEnergyLoss() << std::endl << " ParticleType = " << hit.getParticleType() << std::endl << " Pabs = " << hit.getPabs() << std::endl - << " Energy of primary particle (ID = " << hit.getTrackID() << ") = " << hit.getIncidentEnergy() << " (MeV)" - << std::endl + << " Energy of primary particle (ID = " << hit.getTrackID() << " hit type = " << hit.getProdType() + << ") = " << hit.getIncidentEnergy() << " (MeV)" << std::endl << " Entry point in Bsc unit number " << hit.getUnitID() << " is: " << hit.getEntry() << " (mm)" << std::endl; os << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; return os; diff --git a/SimG4CMS/Forward/src/MtdSD.cc b/SimG4CMS/Forward/src/MtdSD.cc index fc5869c86b363..3bf9728baad54 100644 --- a/SimG4CMS/Forward/src/MtdSD.cc +++ b/SimG4CMS/Forward/src/MtdSD.cc @@ -16,8 +16,6 @@ #include -using namespace MtdHitCategory; - //------------------------------------------------------------------- MtdSD::MtdSD(const std::string& name, const SensitiveDetectorCatalog& clg, @@ -111,32 +109,12 @@ int MtdSD::getTrackID(const G4Track* aTrack) { if (!trkInfo->storeTrack()) { theID = trkInfo->idLastStoredAncestor(); } - if (theID >= static_cast(PSimHit::k_tidOffset)) { - edm::LogError("MtdSim") << " SimTrack ID " << theID << " exceeds maximum allowed by PSimHit identifier" - << PSimHit::k_tidOffset << " unreliable MTD hit type"; - } - if (rname == "FastTimerRegionSensBTL") { - if (trkInfo->isInTrkFromBackscattering()) { - theID = PSimHit::addTrackIdOffset(theID, k_idFromCaloOffset); - } else if (trkInfo->isExtSecondary() && !trkInfo->isInTrkFromBackscattering() && !trkInfo->storeTrack()) { - theID = PSimHit::addTrackIdOffset(theID, k_idsecOffset); - } else if (trkInfo->isBTLlooper()) { - theID = PSimHit::addTrackIdOffset(theID, k_idloopOffset); - } -#ifdef EDM_ML_DEBUG - edm::LogVerbatim("MtdSim") << "MtdSD: Track ID: " << aTrack->GetTrackID() - << " BTL Track ID: " << trkInfo->mcTruthID() << ":" << theID; -#endif - } else if (rname == "FastTimerRegionSensETL") { - if (hitClassID == k_idETLfromBack) { - theID = PSimHit::addTrackIdOffset(theID, k_idETLfromBack); - } #ifdef EDM_ML_DEBUG - edm::LogVerbatim("MtdSim") << "MtdSD: Track ID: " << aTrack->GetTrackID() - << " ETL Track ID: " << trkInfo->mcTruthID() << ":" << theID; + edm::LogVerbatim("MtdSim") << "MtdSD: current Track ID: " << aTrack->GetTrackID() + << " stored Track ID: " << trkInfo->mcTruthID() << ":" << theID; #endif - // In the case of ECAL GFlash fast spot may be inside MTD and should be ignored - } else { + // In the case of ECAL GFlash fast spot may be inside MTD and should be ignored + if (rname != "FastTimerRegionSensBTL" && rname != "FastTimerRegionSensETL") { throw cms::Exception("MtdSDError") << "MtdSD called in incorrect region " << rname; } } else { @@ -148,20 +126,31 @@ int MtdSD::getTrackID(const G4Track* aTrack) { } void MtdSD::setHitClassID(const G4Step* aStep) { + hitClassID = 0; TrackInformation* trkInfo = cmsTrackInformation(aStep->GetTrack()); if (nullptr == trkInfo) { return; } const G4String& rname = aStep->GetTrack()->GetVolume()->GetLogicalVolume()->GetRegion()->GetName(); - if (rname == "FastTimerRegionSensETL") { + if (rname == "FastTimerRegionSensBTL") { + if (trkInfo->isInTrkFromBackscattering()) { + hitClassID = SimHitCategory::prodTypeMTD[3]; + } else if (trkInfo->isExtSecondary() && !trkInfo->isInTrkFromBackscattering() && !trkInfo->storeTrack()) { + hitClassID = SimHitCategory::prodTypeMTD[1]; + } else if (trkInfo->isBTLlooper()) { + hitClassID = SimHitCategory::prodTypeMTD[2]; + } + } else if (rname == "FastTimerRegionSensETL") { double zin = std::abs(aStep->GetPreStepPoint()->GetPosition().z()); double zout = std::abs(aStep->GetPostStepPoint()->GetPosition().z()); if (zout - zin < 0.) { - hitClassID = k_idETLfromBack; + hitClassID = SimHitCategory::prodTypeMTD[4]; trkInfo->setETLfromBack(); } else { - hitClassID = 0; trkInfo->setETLfromFront(); } } +#ifdef EDM_ML_DEBUG + edm::LogVerbatim("MtdSim") << "MtdSD: process type = " << hitClassID; +#endif } diff --git a/SimG4CMS/Forward/src/TimingSD.cc b/SimG4CMS/Forward/src/TimingSD.cc index 39c90a5c7a4ac..3998b657d318c 100644 --- a/SimG4CMS/Forward/src/TimingSD.cc +++ b/SimG4CMS/Forward/src/TimingSD.cc @@ -238,6 +238,7 @@ bool TimingSD::checkHit(const G4Step*, BscG4Hit* hit) { hit->setParentId(theTrack->GetParentID()); hit->setProcessId(theEnumerator->processId(theTrack->GetCreatorProcess())); + hit->setHitProdType(hitClassID); hit->setVertexPosition(theTrack->GetVertexPosition()); } @@ -294,6 +295,7 @@ void TimingSD::createNewHit(const G4Step* aStep) { currentHit->setParentId(theTrack->GetParentID()); currentHit->setProcessId(theEnumerator->processId(theTrack->GetCreatorProcess())); + currentHit->setHitProdType(hitClassID); currentHit->setVertexPosition(theTrack->GetVertexPosition()); diff --git a/SimG4Core/Application/test/SimHitCaloHitDumper.cc b/SimG4Core/Application/test/SimHitCaloHitDumper.cc index de105b3387768..80eeb75b26ba5 100644 --- a/SimG4Core/Application/test/SimHitCaloHitDumper.cc +++ b/SimG4Core/Application/test/SimHitCaloHitDumper.cc @@ -373,8 +373,7 @@ void SimHitCaloHitDumper::analyze(const edm::Event& iEvent, const edm::EventSetu for (int ihit = 0; ihit < (*icoll).first; ++ihit) { edm::LogPrint("SimHitCaloHitDumper") << theMTDHits[nhit] << " Energy = " << theMTDHits[nhit].energyLoss() - << " tid orig/offset= " << theMTDHits[nhit].originalTrackId() << " " << theMTDHits[nhit].offsetTrackId() - << " Track Id = " << theMTDHits[nhit].trackId(); + << " Track Id = " << theMTDHits[nhit].trackId() << " hit type = " << theMTDHits[nhit].hitProdType(); nhit++; } } diff --git a/SimG4Core/Notification/src/SimTrackManager.cc b/SimG4Core/Notification/src/SimTrackManager.cc index 16162236a34ea..73d771f0e6128 100644 --- a/SimG4Core/Notification/src/SimTrackManager.cc +++ b/SimG4Core/Notification/src/SimTrackManager.cc @@ -175,10 +175,6 @@ void SimTrackManager::reallyStoreTracks() { } } - if (id >= static_cast(PSimHit::k_tidOffset)) { - edm::LogWarning("SimTrackManager::reallyStoreTracks") - << " SimTrack ID " << id << " exceeds maximum allowed by PSimHit identifier" << PSimHit::k_tidOffset; - } TmpSimTrack* g4simtrack = new TmpSimTrack(id, trkH->particleID(), trkH->momentum(), trkH->totalEnergy(), ivertex, ig, pm, spos, smom); g4simtrack->copyCrossedBoundaryVars(trkH); diff --git a/SimGeneral/CaloAnalysis/plugins/MtdTruthAccumulator.cc b/SimGeneral/CaloAnalysis/plugins/MtdTruthAccumulator.cc index 930a1b5daec5e..0def9d5d6f417 100644 --- a/SimGeneral/CaloAnalysis/plugins/MtdTruthAccumulator.cc +++ b/SimGeneral/CaloAnalysis/plugins/MtdTruthAccumulator.cc @@ -42,6 +42,7 @@ #include "SimDataFormats/CaloAnalysis/interface/MtdSimTracksterFwd.h" #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h" #include "SimDataFormats/TrackingHit/interface/PSimHit.h" +#include "SimDataFormats/TrackingHit/interface/SimHitCategory.h" #include "SimDataFormats/Vertex/interface/SimVertex.h" #include "SimGeneral/MixingModule/interface/DecayGraph.h" @@ -59,14 +60,17 @@ #include "Geometry/MTDGeometryBuilder/interface/ProxyMTDTopology.h" #include "Geometry/MTDGeometryBuilder/interface/RectangularMTDTopology.h" -#include "SimG4CMS/Forward/interface/MtdHitCategory.h" - #include namespace { using Index_t = unsigned; - using Barcode_t = int; + using Barcode_t = int64_t; const std::string messageCategoryGraph_("MtdTruthAccumulatorGraphProducer"); + inline int64_t tidAndProc(int a, unsigned short b) { + auto tmpout = static_cast(a) << 32; + tmpout |= b; + return tmpout; + } } // namespace class MtdTruthAccumulator : public DigiAccumulatorMixMod { @@ -88,10 +92,11 @@ class MtdTruthAccumulator : public DigiAccumulatorMixMod { /** @brief Fills the supplied vector with pointers to the SimHits, checking * for bad modules if required */ template - void fillSimHits(std::vector> &returnValue, - std::unordered_map>> &simTrackDetIdMap, - const T &event, - const edm::EventSetup &setup); + void fillSimHits( + std::vector> &returnValue, + std::unordered_map>> &simTrackDetIdMap, + const T &event, + const edm::EventSetup &setup); const std::string messageCategory_; @@ -164,13 +169,14 @@ class MtdTruthAccumulator : public DigiAccumulatorMixMod { /* Graph utility functions */ namespace { + class CaloParticle_dfs_visitor : public boost::default_dfs_visitor { public: CaloParticle_dfs_visitor( MtdTruthAccumulator::OutputCollections &output, MtdTruthAccumulator::calo_particles &caloParticles, std::unordered_multimap &simHitBarcodeToIndex, - std::unordered_map>> &simTrackDetIdMap, + std::unordered_map>> &simTrackDetIdMap, std::unordered_map &vertex_time_map, Selector selector) : output_(output), @@ -188,10 +194,9 @@ namespace { auto const vertex_property = get(vertex_name, g, u); if (!vertex_property.simTrack) return; - // -- loop over possible trackIdOffsets to save also sim clusters from non-direct hits - for (unsigned int offset = 0; offset < MtdHitCategory::n_categories + 1; offset++) { - auto trackIdx = vertex_property.simTrack->trackId(); - trackIdx += offset * (static_cast(PSimHit::k_tidOffset)); + // -- loop over possible hitProdTypes to save also sim clusters from non-direct hits + for (unsigned int offset = 0; offset < SimHitCategory::nCategoriesMTD; offset++) { + auto trackIdx = tidAndProc(vertex_property.simTrack->trackId(), offset); IfLogDebug(DEBUG, messageCategoryGraph_) << " Found " << simHitBarcodeToIndex_.count(trackIdx) << " associated simHits" << std::endl; if (simHitBarcodeToIndex_.count(trackIdx)) { @@ -205,12 +210,10 @@ namespace { simcluster.addHitAndFraction(hit_and_energy.first, hit_and_energy.second); simcluster.addHitEnergy(hit_and_energy.second); simcluster.addHitTime(std::get<1>( - simTrackDetIdMap_[simcluster.g4Tracks()[0].trackId() + - offset * (static_cast(PSimHit::k_tidOffset))][hit_and_energy.first])); + simTrackDetIdMap_[tidAndProc(simcluster.g4Tracks()[0].trackId(), offset)][hit_and_energy.first])); simcluster.addHitPosition(std::get<2>( - simTrackDetIdMap_[simcluster.g4Tracks()[0].trackId() + - offset * (static_cast(PSimHit::k_tidOffset))][hit_and_energy.first])); - simcluster.setTrackIdOffset(offset); + simTrackDetIdMap_[tidAndProc(simcluster.g4Tracks()[0].trackId(), offset)][hit_and_energy.first])); + simcluster.setHitProdType(offset); } } } @@ -248,7 +251,7 @@ namespace { MtdTruthAccumulator::OutputCollections &output_; MtdTruthAccumulator::calo_particles &caloParticles_; std::unordered_multimap &simHitBarcodeToIndex_; - std::unordered_map>> &simTrackDetIdMap_; + std::unordered_map>> &simTrackDetIdMap_; std::unordered_map &vertex_time_map_; Selector selector_; }; @@ -472,7 +475,7 @@ void MtdTruthAccumulator::finalizeEvent(edm::Event &event, edm::EventSetup const SimLCz = 0.; tmpLC.addCluIndex(SC_index); tmpLC.computeClusterTime(); - tmpLC.setTrackIdOffset(sc.trackIdOffset()); // add trackIdoffset + tmpLC.setHitProdType(sc.hitProdType()); // add hitProdType output_.pMtdSimLayerClusters->push_back(tmpLC); LC_indices.push_back(LC_index); LC_index++; @@ -605,13 +608,14 @@ void MtdTruthAccumulator::accumulateEvent(const T &event, event.getByLabel(genParticleLabel_, hGenParticleIndices); std::vector> simHitPointers; - std::unordered_map>> simTrackDetIdMap; + std::unordered_map>> simTrackDetIdMap; fillSimHits(simHitPointers, simTrackDetIdMap, event, setup); // Clear maps from previous event fill them for this one m_simHitBarcodeToIndex.clear(); for (unsigned int i = 0; i < simHitPointers.size(); ++i) { - m_simHitBarcodeToIndex.emplace(simHitPointers[i].second->trackId(), i); + m_simHitBarcodeToIndex.emplace( + tidAndProc(simHitPointers[i].second->trackId(), simHitPointers[i].second->hitProdType()), i); } auto const &tracks = *hSimTracks; @@ -740,7 +744,7 @@ void MtdTruthAccumulator::accumulateEvent(const T &event, template void MtdTruthAccumulator::fillSimHits( std::vector> &returnValue, - std::unordered_map>> &simTrackDetIdMap, + std::unordered_map>> &simTrackDetIdMap, const T &event, const edm::EventSetup &setup) { using namespace geant_units::operators; @@ -781,19 +785,20 @@ void MtdTruthAccumulator::fillSimHits( uniqueId |= pixel.first << 16; uniqueId |= pixel.second; - std::get<0>(simTrackDetIdMap[simHit.trackId()][uniqueId]) += simHit.energyLoss(); + int64_t trackIdx = tidAndProc(simHit.trackId(), simHit.hitProdType()); + std::get<0>(simTrackDetIdMap[trackIdx][uniqueId]) += simHit.energyLoss(); m_detIdToTotalSimEnergy[uniqueId] += simHit.energyLoss(); // --- Get the time of the first SIM hit in the cell - if (std::get<1>(simTrackDetIdMap[simHit.trackId()][uniqueId]) == 0. || - simHit.tof() < std::get<1>(simTrackDetIdMap[simHit.trackId()][uniqueId])) { - std::get<1>(simTrackDetIdMap[simHit.trackId()][uniqueId]) = simHit.tof(); + if (std::get<1>(simTrackDetIdMap[trackIdx][uniqueId]) == 0. || + simHit.tof() < std::get<1>(simTrackDetIdMap[trackIdx][uniqueId])) { + std::get<1>(simTrackDetIdMap[trackIdx][uniqueId]) = simHit.tof(); } - float xSim = std::get<2>(simTrackDetIdMap[simHit.trackId()][uniqueId]).x() + simscaled.x() * simHit.energyLoss(); - float ySim = std::get<2>(simTrackDetIdMap[simHit.trackId()][uniqueId]).y() + simscaled.y() * simHit.energyLoss(); - float zSim = std::get<2>(simTrackDetIdMap[simHit.trackId()][uniqueId]).z() + simscaled.z() * simHit.energyLoss(); + float xSim = std::get<2>(simTrackDetIdMap[trackIdx][uniqueId]).x() + simscaled.x() * simHit.energyLoss(); + float ySim = std::get<2>(simTrackDetIdMap[trackIdx][uniqueId]).y() + simscaled.y() * simHit.energyLoss(); + float zSim = std::get<2>(simTrackDetIdMap[trackIdx][uniqueId]).z() + simscaled.z() * simHit.energyLoss(); LocalPoint posSim(xSim, ySim, zSim); - std::get<2>(simTrackDetIdMap[simHit.trackId()][uniqueId]) = posSim; + std::get<2>(simTrackDetIdMap[trackIdx][uniqueId]) = posSim; #ifdef PRINT_DEBUG IfLogDebug(DEBUG, messageCategory_) diff --git a/Validation/MtdValidation/plugins/BtlLocalRecoValidation.cc b/Validation/MtdValidation/plugins/BtlLocalRecoValidation.cc index 219a6fd5532a2..1beab1d520715 100644 --- a/Validation/MtdValidation/plugins/BtlLocalRecoValidation.cc +++ b/Validation/MtdValidation/plugins/BtlLocalRecoValidation.cc @@ -686,7 +686,7 @@ void BtlLocalRecoValidation::analyze(const edm::Event& iEvent, const edm::EventS float simClusTime = (*simClusterRef).simLCTime(); LocalPoint simClusLocalPos = (*simClusterRef).simLCPos(); const auto& simClusGlobalPos = genericDet->toGlobal(simClusLocalPos); - unsigned int idOffset = (*simClusterRef).trackIdOffset(); + unsigned int idOffset = (*simClusterRef).hitProdType(); float time_res = cluster.time() - simClusTime; float energy_res = cluster.energy() - simClusEnergy; diff --git a/Validation/MtdValidation/plugins/BtlSimHitsValidation.cc b/Validation/MtdValidation/plugins/BtlSimHitsValidation.cc index edff176f70e3c..f4741ca1893dc 100644 --- a/Validation/MtdValidation/plugins/BtlSimHitsValidation.cc +++ b/Validation/MtdValidation/plugins/BtlSimHitsValidation.cc @@ -41,6 +41,14 @@ #include "MTDHit.h" +namespace { + + inline uint64_t cantor(uint64_t a, uint64_t b) { return (a + b + 1) * (a + b) / 2 + b; } + + inline uint64_t hash3(uint64_t a, uint64_t b, uint64_t c) { return cantor(a, cantor(b, c)); } + +} // namespace + class BtlSimHitsValidation : public DQMEDAnalyzer { public: explicit BtlSimHitsValidation(const edm::ParameterSet&); @@ -181,26 +189,28 @@ void BtlSimHitsValidation::analyze(const edm::Event& iEvent, const edm::EventSet DetId id = simHit.detUnitId(); // --- Build a global track ID by combining the SIM hit's eventId and trackId - uint64_t globalTrkID = ((uint64_t)simHit.eventId().rawId() << 32) | simHit.trackId(); + uint64_t cmbID = hash3(simHit.eventId().rawId(), simHit.trackId(), simHit.hitProdType()); // --- Define geoId from MTDTopology to group cells with hits by sensor module BTLDetId detId(id.rawId()); DetId geoId = detId.geographicalId(MTDTopologyMode::crysLayoutFromTopoMode(topology->getMTDTopologyMode())); // --- Sum the energies of SIM hits with the same track ID in the same cell - m_btlHitsPerCellAndModule[geoId.rawId()][id.rawId()][globalTrkID].energy += convertGeVToMeV(simHit.energyLoss()); + m_btlHitsPerCellAndModule[geoId.rawId()][id.rawId()][cmbID].energy += convertGeVToMeV(simHit.energyLoss()); // --- Assign the time and position of the first SIM hit in time to the accumulated hit - if (m_btlHitsPerCellAndModule[geoId.rawId()][id.rawId()][globalTrkID].time == 0. || - simHit.tof() < m_btlHitsPerCellAndModule[geoId.rawId()][id.rawId()][globalTrkID].time) { - m_btlHitsPerCellAndModule[geoId.rawId()][id.rawId()][globalTrkID].time = simHit.tof(); + if (m_btlHitsPerCellAndModule[geoId.rawId()][id.rawId()][cmbID].time == 0. || + simHit.tof() < m_btlHitsPerCellAndModule[geoId.rawId()][id.rawId()][cmbID].time) { + m_btlHitsPerCellAndModule[geoId.rawId()][id.rawId()][cmbID].time = simHit.tof(); auto hit_pos = simHit.localPosition(); - m_btlHitsPerCellAndModule[geoId.rawId()][id.rawId()][globalTrkID].x = hit_pos.x(); - m_btlHitsPerCellAndModule[geoId.rawId()][id.rawId()][globalTrkID].y = hit_pos.y(); - m_btlHitsPerCellAndModule[geoId.rawId()][id.rawId()][globalTrkID].z = hit_pos.z(); + m_btlHitsPerCellAndModule[geoId.rawId()][id.rawId()][cmbID].x = hit_pos.x(); + m_btlHitsPerCellAndModule[geoId.rawId()][id.rawId()][cmbID].y = hit_pos.y(); + m_btlHitsPerCellAndModule[geoId.rawId()][id.rawId()][cmbID].z = hit_pos.z(); - m_btlHitsPerCellAndModule[geoId.rawId()][id.rawId()][globalTrkID].thetaAtEntry = simHit.thetaAtEntry(); + m_btlHitsPerCellAndModule[geoId.rawId()][id.rawId()][cmbID].thetaAtEntry = simHit.thetaAtEntry(); + m_btlHitsPerCellAndModule[geoId.rawId()][id.rawId()][cmbID].prodType = + 2 * simHit.hitProdType(); // for backward compatibility multiply x 2 } } // simHit loop @@ -217,6 +227,7 @@ void BtlSimHitsValidation::analyze(const edm::Event& iEvent, const edm::EventSet meNhits_->Fill(log10(Nhits)); // --- Loop over the BTL sensor modules + std::vector v_hitID; for (const auto& module : m_btlHitsPerCellAndModule) { // --- Loop over the BTL cells for (auto const& cell : module.second) { @@ -229,7 +240,7 @@ void BtlSimHitsValidation::analyze(const edm::Event& iEvent, const edm::EventSet } // --- Loop over the hits in the cell, sum the hit energies and store the hit IDs in a vector - std::vector v_hitID; + v_hitID.clear(); float ene_tot_cell = 0.; for (auto const& hit : m_hits) { @@ -319,8 +330,8 @@ void BtlSimHitsValidation::analyze(const edm::Event& iEvent, const edm::EventSet meNtrkPerCell_->Fill(v_hitID.size()); // --- First hit in the cell - int trackID = (int)(v_hitID[0] & 0xFFFFFFFF) / 100000000; - meHitTrkID1_->Fill(trackID); + int prodID = (int)(m_hits.at(v_hitID[0]).prodType); + meHitTrkID1_->Fill(prodID); // cells in the bulk of energy distribution if (ene_tot_cell < cellEneCut_) { @@ -333,8 +344,8 @@ void BtlSimHitsValidation::analyze(const edm::Event& iEvent, const edm::EventSet // --- Second hit in the cell if (v_hitID.size() == 2) { - trackID = (int)(v_hitID[1] & 0xFFFFFFFF) / 100000000; - meHitTrkID2_->Fill(trackID); + prodID = (int)(m_hits.at(v_hitID[1]).prodType); + meHitTrkID2_->Fill(prodID); meHitDeltaT2_->Fill(deltaT_max); @@ -353,8 +364,8 @@ void BtlSimHitsValidation::analyze(const edm::Event& iEvent, const edm::EventSet } // --- Third hit in the cell else if (v_hitID.size() == 3) { - trackID = (int)(v_hitID[2] & 0xFFFFFFFF) / 100000000; - meHitTrkID3_->Fill(trackID); + prodID = (int)(m_hits.at(v_hitID[2]).prodType); + meHitTrkID3_->Fill(prodID); meHitDeltaT3_->Fill(deltaT_max); @@ -371,8 +382,8 @@ void BtlSimHitsValidation::analyze(const edm::Event& iEvent, const edm::EventSet // --- Fourth hit in the cell and next ones else if (v_hitID.size() >= 4) { for (unsigned int ihit = 3; ihit < v_hitID.size(); ++ihit) { - trackID = (int)(v_hitID[ihit] & 0xFFFFFFFF) / 100000000; - meHitTrkID4_->Fill(trackID); + prodID = (int)(m_hits.at(v_hitID[ihit]).prodType); + meHitTrkID4_->Fill(prodID); // cells in the bulk of energy distribution if (ene_tot_cell < cellEneCut_) { diff --git a/Validation/MtdValidation/plugins/EtlLocalRecoValidation.cc b/Validation/MtdValidation/plugins/EtlLocalRecoValidation.cc index 42b95272511d8..f0a36f541d639 100644 --- a/Validation/MtdValidation/plugins/EtlLocalRecoValidation.cc +++ b/Validation/MtdValidation/plugins/EtlLocalRecoValidation.cc @@ -526,7 +526,7 @@ void EtlLocalRecoValidation::analyze(const edm::Event& iEvent, const edm::EventS (*itp.first).second; // the range of itp.first, itp.second should be always 1 for (unsigned int i = 0; i < simClustersRefs.size(); i++) { const auto& simClusterRef = simClustersRefs[i]; - unsigned int idOffset = (*simClusterRef).trackIdOffset(); + unsigned int idOffset = (*simClusterRef).hitProdType(); meCluTrackIdOffset_[idet]->Fill(float(idOffset)); diff --git a/Validation/MtdValidation/plugins/EtlSimHitsValidation.cc b/Validation/MtdValidation/plugins/EtlSimHitsValidation.cc index 69e926b7cb5e1..90f3e38341698 100644 --- a/Validation/MtdValidation/plugins/EtlSimHitsValidation.cc +++ b/Validation/MtdValidation/plugins/EtlSimHitsValidation.cc @@ -182,7 +182,7 @@ void EtlSimHitsValidation::analyze(const edm::Event& iEvent, const edm::EventSet (simHitIt->second).y = hit_pos.y(); (simHitIt->second).z = hit_pos.z(); - if (simHit.offsetTrackId() == 0) { + if (simHit.hitProdType() == 0) { if (simHit.exitPoint() != simHit.entryPoint()) { (simHitIt->second).thetaAtEntry = angle_units::operators::convertRadToDeg((simHit.exitPoint() - simHit.entryPoint()).bareTheta()); diff --git a/Validation/MtdValidation/plugins/MTDHit.h b/Validation/MtdValidation/plugins/MTDHit.h index c3186d96d8b93..aa67ed716e83e 100644 --- a/Validation/MtdValidation/plugins/MTDHit.h +++ b/Validation/MtdValidation/plugins/MTDHit.h @@ -8,6 +8,7 @@ struct MTDHit { float y; float z; float thetaAtEntry; + unsigned short prodType; }; #endif //Validation_MtdValidation_MTDHit_h diff --git a/Validation/MtdValidation/plugins/MtdTracksValidation.cc b/Validation/MtdValidation/plugins/MtdTracksValidation.cc index a8f966f89c869..6e8789bcebb75 100644 --- a/Validation/MtdValidation/plugins/MtdTracksValidation.cc +++ b/Validation/MtdValidation/plugins/MtdTracksValidation.cc @@ -906,13 +906,13 @@ void MtdTracksValidation::analyze(const edm::Event& iEvent, const edm::EventSetu // Find the first direct hit in time directSimClusIt = std::find_if(simClustersRefs.begin(), simClustersRefs.end(), [](const auto& simCluster) { MTDDetId mtddetid = simCluster->detIds_and_rows().front().first; - return (mtddetid.mtdSubDetector() == 1 && simCluster->trackIdOffset() == 0); + return (mtddetid.mtdSubDetector() == 1 && simCluster->hitProdType() == 0); }); // Check if TP has direct or other sim cluster for BTL for (const auto& simClusterRef : simClustersRefs) { if (directSimClusIt != simClustersRefs.end() && simClusterRef == *directSimClusIt) { isTPmtdDirectBTL = true; - } else if (simClusterRef->trackIdOffset() != 0) { + } else if (simClusterRef->hitProdType() != 0) { isTPmtdOtherBTL = true; } } @@ -950,7 +950,7 @@ void MtdTracksValidation::analyze(const edm::Event& iEvent, const edm::EventSetu GlobalPoint simClusRecoMatchGlobalPos = geomUtil.globalPosition(detidRecoMatch, simClusRecoMatchLocalPos); - simClusterRef_RecoMatch_trackIdOff = simClusterRef_RecoMatch->trackIdOffset(); + simClusterRef_RecoMatch_trackIdOff = simClusterRef_RecoMatch->hitProdType(); simClusterRef_RecoMatch_DeltaZ = simClusRecoMatchGlobalPos.z() - simClusGlobalPos.z(); simClusterRef_RecoMatch_DeltaPhi = simClusRecoMatchGlobalPos.phi() - simClusGlobalPos.phi(); simClusterRef_RecoMatch_DeltaT = simClusterRef_RecoMatch->simLCTime() - directSimClus->simLCTime(); @@ -964,7 +964,7 @@ void MtdTracksValidation::analyze(const edm::Event& iEvent, const edm::EventSetu isTPmtdDirectCorrectBTL = true; simClusterEarliestTime_correctAssoc = earliestSimHitTimeInSimCluster((*simClusterIt)); - } else if (simClusterRef_RecoMatch->trackIdOffset() != 0) { + } else if (simClusterRef_RecoMatch->hitProdType() != 0) { isTPmtdOtherCorrectBTL = true; } } @@ -1374,7 +1374,7 @@ void MtdTracksValidation::analyze(const edm::Event& iEvent, const edm::EventSetu meBTLTrackMatchedTPnomtdAssocTrackID_->Fill(1); } } - meBTLTrackMatchedTPnomtdAssocTrackIdOff_->Fill(sc->trackIdOffset()); + meBTLTrackMatchedTPnomtdAssocTrackIdOff_->Fill(sc->hitProdType()); } } }