Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-05-29 03:17:58 

0001 // Author: Felice Pantaleo, felice.pantaleo@cern.ch 08/2024
0002 
0003 #include "FWCore/Framework/interface/global/EDProducer.h"
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/MakerMacros.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "FWCore/Utilities/interface/EDGetToken.h"
0009 #include "DataFormats/HGCalReco/interface/Trackster.h"
0010 #include "SimDataFormats/Associations/interface/TICLAssociationMap.h"
0011 #include "DataFormats/Provenance/interface/ProductID.h"
0012 #include "FWCore/Framework/interface/ESHandle.h"
0013 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
0014 
0015 class AllTracksterToSimTracksterAssociatorsByLCsProducer : public edm::global::EDProducer<> {
0016 public:
0017   explicit AllTracksterToSimTracksterAssociatorsByLCsProducer(const edm::ParameterSet&);
0018   ~AllTracksterToSimTracksterAssociatorsByLCsProducer() override = default;
0019 
0020   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0021 
0022 private:
0023   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0024 
0025   std::vector<std::pair<std::string, edm::EDGetTokenT<std::vector<ticl::Trackster>>>> tracksterCollectionTokens_;
0026   std::vector<std::pair<std::string, edm::EDGetTokenT<std::vector<ticl::Trackster>>>> simTracksterCollectionTokens_;
0027   edm::EDGetTokenT<std::vector<reco::CaloCluster>> layerClustersToken_;
0028   std::vector<std::pair<
0029       std::string,
0030       edm::EDGetTokenT<
0031           ticl::AssociationMap<ticl::mapWithSharedEnergy, std::vector<reco::CaloCluster>, std::vector<ticl::Trackster>>>>>
0032       layerClusterToTracksterMapTokens_;
0033   std::vector<std::pair<
0034       std::string,
0035       edm::EDGetTokenT<
0036           ticl::AssociationMap<ticl::mapWithSharedEnergy, std::vector<reco::CaloCluster>, std::vector<ticl::Trackster>>>>>
0037       layerClusterToSimTracksterMapTokens_;
0038 };
0039 
0040 AllTracksterToSimTracksterAssociatorsByLCsProducer::AllTracksterToSimTracksterAssociatorsByLCsProducer(
0041     const edm::ParameterSet& pset)
0042     : layerClustersToken_(consumes<std::vector<reco::CaloCluster>>(pset.getParameter<edm::InputTag>("layerClusters"))) {
0043   const auto& tracksterCollections = pset.getParameter<std::vector<edm::InputTag>>("tracksterCollections");
0044 
0045   std::string allLCtoTSAccoc = pset.getParameter<std::string>("allLCtoTSAccoc");
0046   for (const auto& tag : tracksterCollections) {
0047     std::string label = tag.label();
0048     if (!tag.instance().empty()) {
0049       label += tag.instance();
0050     }
0051     tracksterCollectionTokens_.emplace_back(label, consumes<std::vector<ticl::Trackster>>(tag));
0052     layerClusterToTracksterMapTokens_.emplace_back(
0053         label,
0054         consumes<ticl::AssociationMap<ticl::mapWithSharedEnergy,
0055                                       std::vector<reco::CaloCluster>,
0056                                       std::vector<ticl::Trackster>>>(edm::InputTag(allLCtoTSAccoc, label)));
0057   }
0058 
0059   const auto& simTracksterCollections = pset.getParameter<std::vector<edm::InputTag>>("simTracksterCollections");
0060   for (const auto& tag : simTracksterCollections) {
0061     std::string label = tag.label();
0062     if (!tag.instance().empty()) {
0063       label += tag.instance();
0064     }
0065     simTracksterCollectionTokens_.emplace_back(label, consumes<std::vector<ticl::Trackster>>(tag));
0066     layerClusterToSimTracksterMapTokens_.emplace_back(
0067         label,
0068         consumes<ticl::AssociationMap<ticl::mapWithSharedEnergy,
0069                                       std::vector<reco::CaloCluster>,
0070                                       std::vector<ticl::Trackster>>>(edm::InputTag(allLCtoTSAccoc, label)));
0071   }
0072 
0073   // Produce separate association maps for each trackster-simTrackster combination
0074   for (const auto& tracksterToken : tracksterCollectionTokens_) {
0075     for (const auto& simTracksterToken : simTracksterCollectionTokens_) {
0076       std::string instanceLabel = tracksterToken.first + "To" + simTracksterToken.first;
0077       produces<ticl::AssociationMap<ticl::mapWithSharedEnergyAndScore,
0078                                     std::vector<ticl::Trackster>,
0079                                     std::vector<ticl::Trackster>>>(instanceLabel);
0080       std::string reverseInstanceLabel = simTracksterToken.first + "To" + tracksterToken.first;
0081       produces<ticl::AssociationMap<ticl::mapWithSharedEnergyAndScore,
0082                                     std::vector<ticl::Trackster>,
0083                                     std::vector<ticl::Trackster>>>(reverseInstanceLabel);
0084     }
0085   }
0086 }
0087 
0088 void AllTracksterToSimTracksterAssociatorsByLCsProducer::produce(edm::StreamID,
0089                                                                  edm::Event& iEvent,
0090                                                                  const edm::EventSetup&) const {
0091   using namespace edm;
0092 
0093   // Retrieve layer clusters with protection
0094   const auto& layerClustersHandle = iEvent.getHandle(layerClustersToken_);
0095 
0096   // If layer clusters are missing, produce empty maps and return
0097   if (!layerClustersHandle.isValid()) {
0098     edm::LogWarning("MissingInput") << "Layer clusters collection not found. Producing empty maps.";
0099     for (const auto& tracksterToken : tracksterCollectionTokens_) {
0100       for (const auto& simTracksterToken : simTracksterCollectionTokens_) {
0101         iEvent.put(std::make_unique<ticl::AssociationMap<ticl::mapWithSharedEnergyAndScore,
0102                                                          std::vector<ticl::Trackster>,
0103                                                          std::vector<ticl::Trackster>>>(),
0104                    tracksterToken.first + "To" + simTracksterToken.first);
0105 
0106         iEvent.put(std::make_unique<ticl::AssociationMap<ticl::mapWithSharedEnergyAndScore,
0107                                                          std::vector<ticl::Trackster>,
0108                                                          std::vector<ticl::Trackster>>>(),
0109                    simTracksterToken.first + "To" + tracksterToken.first);
0110       }
0111     }
0112     return;
0113   }
0114 
0115   const auto& layerClusters = *layerClustersHandle;
0116 
0117   for (const auto& tracksterToken : tracksterCollectionTokens_) {
0118     Handle<std::vector<ticl::Trackster>> recoTrackstersHandle;
0119     iEvent.getByToken(tracksterToken.second, recoTrackstersHandle);
0120 
0121     if (!recoTrackstersHandle.isValid()) {
0122       edm::LogWarning("MissingInput") << "trackster  collection not found. Producing empty maps.";
0123       for (const auto& simTracksterToken : simTracksterCollectionTokens_) {
0124         iEvent.put(std::make_unique<ticl::AssociationMap<ticl::mapWithSharedEnergyAndScore,
0125                                                          std::vector<ticl::Trackster>,
0126                                                          std::vector<ticl::Trackster>>>(),
0127                    tracksterToken.first + "To" + simTracksterToken.first);
0128 
0129         iEvent.put(std::make_unique<ticl::AssociationMap<ticl::mapWithSharedEnergyAndScore,
0130                                                          std::vector<ticl::Trackster>,
0131                                                          std::vector<ticl::Trackster>>>(),
0132                    simTracksterToken.first + "To" + tracksterToken.first);
0133       }
0134       continue;
0135     }
0136 
0137     const auto& recoTracksters = *recoTrackstersHandle;
0138 
0139     // Retrieve the correct LayerClusterToTracksterMap for the current trackster collection
0140     Handle<ticl::AssociationMap<ticl::mapWithSharedEnergy, std::vector<reco::CaloCluster>, std::vector<ticl::Trackster>>>
0141         layerClusterToTracksterMapHandle;
0142     auto tracksterMapTokenIter =
0143         std::find_if(layerClusterToTracksterMapTokens_.begin(),
0144                      layerClusterToTracksterMapTokens_.end(),
0145                      [&tracksterToken](const auto& pair) { return pair.first == tracksterToken.first; });
0146     if (tracksterMapTokenIter != layerClusterToTracksterMapTokens_.end()) {
0147       iEvent.getByToken(tracksterMapTokenIter->second, layerClusterToTracksterMapHandle);
0148     }
0149     const auto& layerClusterToTracksterMap = *layerClusterToTracksterMapHandle;
0150 
0151     for (const auto& simTracksterToken : simTracksterCollectionTokens_) {
0152       Handle<std::vector<ticl::Trackster>> simTrackstersHandle;
0153       iEvent.getByToken(simTracksterToken.second, simTrackstersHandle);
0154       const auto& simTracksters = *simTrackstersHandle;
0155 
0156       // Retrieve the correct LayerClusterToSimTracksterMap for the current simTrackster collection
0157       Handle<
0158           ticl::AssociationMap<ticl::mapWithSharedEnergy, std::vector<reco::CaloCluster>, std::vector<ticl::Trackster>>>
0159           layerClusterToSimTracksterMapHandle;
0160       auto simTracksterMapTokenIter =
0161           std::find_if(layerClusterToSimTracksterMapTokens_.begin(),
0162                        layerClusterToSimTracksterMapTokens_.end(),
0163                        [&simTracksterToken](const auto& pair) { return pair.first == simTracksterToken.first; });
0164       if (simTracksterMapTokenIter != layerClusterToSimTracksterMapTokens_.end()) {
0165         iEvent.getByToken(simTracksterMapTokenIter->second, layerClusterToSimTracksterMapHandle);
0166       }
0167       const auto& layerClusterToSimTracksterMap = *layerClusterToSimTracksterMapHandle;
0168 
0169       // Create the association maps
0170       auto tracksterToSimTracksterMap = std::make_unique<ticl::AssociationMap<ticl::mapWithSharedEnergyAndScore,
0171                                                                               std::vector<ticl::Trackster>,
0172                                                                               std::vector<ticl::Trackster>>>(
0173           recoTrackstersHandle, simTrackstersHandle, iEvent);
0174       auto simTracksterToTracksterMap = std::make_unique<ticl::AssociationMap<ticl::mapWithSharedEnergyAndScore,
0175                                                                               std::vector<ticl::Trackster>,
0176                                                                               std::vector<ticl::Trackster>>>(
0177           simTrackstersHandle, recoTrackstersHandle, iEvent);
0178 
0179       for (unsigned int tracksterIndex = 0; tracksterIndex < recoTracksters.size(); ++tracksterIndex) {
0180         const auto& recoTrackster = recoTracksters[tracksterIndex];
0181         edm::Ref<std::vector<ticl::Trackster>> recoTracksterRef(recoTrackstersHandle, tracksterIndex);
0182         const auto& layerClustersIds = recoTrackster.vertices();
0183         float recoToSimScoresDenominator = 0.f;
0184         ticl::AssociationMap<ticl::mapWithSharedEnergy> layerClusterToAssociatedSimTracksterMap(
0185             layerClustersIds.size());
0186         std::vector<unsigned int> associatedSimTracksterIndices;
0187         for (unsigned int i = 0; i < layerClustersIds.size(); ++i) {
0188           unsigned int layerClusterId = layerClustersIds[i];
0189           const auto& layerCluster = layerClusters[layerClusterId];
0190           float recoFraction = 1.f / recoTrackster.vertex_multiplicity(i);
0191           float squaredRecoFraction = recoFraction * recoFraction;
0192           float squaredLayerClusterEnergy = layerCluster.energy() * layerCluster.energy();
0193           recoToSimScoresDenominator += squaredLayerClusterEnergy * squaredRecoFraction;
0194           const auto& simTracksterVec = layerClusterToSimTracksterMap[layerClusterId];
0195           for (const auto& simTracksterElement : simTracksterVec) {
0196             auto simTracksterIndex = simTracksterElement.index();
0197             auto simSharedEnergy = simTracksterElement.sharedEnergy();
0198             layerClusterToAssociatedSimTracksterMap[i].emplace_back(simTracksterIndex, simSharedEnergy);
0199             associatedSimTracksterIndices.push_back(simTracksterIndex);
0200           }
0201         }
0202 
0203         // Keep only unique associatedSimTracksterIndices
0204         std::sort(associatedSimTracksterIndices.begin(), associatedSimTracksterIndices.end());
0205         associatedSimTracksterIndices.erase(
0206             std::unique(associatedSimTracksterIndices.begin(), associatedSimTracksterIndices.end()),
0207             associatedSimTracksterIndices.end());
0208 
0209         // Add missing sim tracksters with 0 shared energy to layerClusterToAssociatedSimTracksterMap
0210         for (unsigned int i = 0; i < layerClustersIds.size(); ++i) {
0211           unsigned int layerClusterId = layerClustersIds[i];
0212           const auto& simTracksterVec = layerClusterToSimTracksterMap[layerClusterId];
0213           for (unsigned int simTracksterIndex : associatedSimTracksterIndices) {
0214             if (std::find_if(simTracksterVec.begin(), simTracksterVec.end(), [simTracksterIndex](const auto& pair) {
0215                   return pair.index() == simTracksterIndex;
0216                 }) == simTracksterVec.end()) {
0217               layerClusterToAssociatedSimTracksterMap[i].push_back({simTracksterIndex, 0.f});
0218             }
0219           }
0220         }
0221 
0222         const float invDenominator = 1.f / recoToSimScoresDenominator;
0223 
0224         for (unsigned int i = 0; i < layerClustersIds.size(); ++i) {
0225           unsigned int layerClusterId = layerClustersIds[i];
0226           const auto& layerCluster = layerClusters[layerClusterId];
0227           float recoFraction = 1.f / recoTrackster.vertex_multiplicity(i);
0228           float squaredLayerClusterEnergy = layerCluster.energy() * layerCluster.energy();
0229           float recoSharedEnergy = layerCluster.energy() * recoFraction;
0230           float invLayerClusterEnergy = 1.f / layerCluster.energy();
0231           const auto& simTracksterVec = layerClusterToAssociatedSimTracksterMap[i];
0232           for (const auto& simTracksterElement : simTracksterVec) {
0233             auto simTracksterIndex = simTracksterElement.index();
0234             float simSharedEnergy = simTracksterElement.sharedEnergy();
0235             edm::Ref<std::vector<ticl::Trackster>> simTracksterRef(simTrackstersHandle, simTracksterIndex);
0236             float sharedEnergy = std::min(simSharedEnergy, recoSharedEnergy);
0237             float simFraction = simSharedEnergy * invLayerClusterEnergy;
0238             /* RecoToSim score logic:
0239              - simFraction >= 0 && recoFraction == 0 : simtrackster contains non-reco associated elements : ignore in recoToSim association
0240              - simFraction == 0 && recoFraction > 0 : rechits not present in sim trackster : penalty in score by recoFraction*E
0241              - simFraction == 1 && recoFraction == 1 : good association
0242              - 1 > recoFraction > simFraction > 0 :  sim does not contain some reco energy, penalty in score by the additional part : (recoFraction-simFraction)*E
0243              - 1 > simFraction> recoFraction > 0 : consider as good association, as there is enough sim to cover the reco
0244             */
0245             float recoMinusSimFraction = std::max(0.f, recoFraction - simFraction);
0246             float score = invDenominator * recoMinusSimFraction * recoMinusSimFraction * squaredLayerClusterEnergy;
0247             tracksterToSimTracksterMap->insert(recoTracksterRef, simTracksterRef, sharedEnergy, score);
0248           }
0249         }
0250       }
0251 
0252       for (unsigned int tracksterIndex = 0; tracksterIndex < simTracksters.size(); ++tracksterIndex) {
0253         const auto& simTrackster = simTracksters[tracksterIndex];
0254         edm::Ref<std::vector<ticl::Trackster>> simTracksterRef(simTrackstersHandle, tracksterIndex);
0255         const auto& layerClustersIds = simTrackster.vertices();
0256         float simToRecoScoresDenominator = 0.f;
0257         ticl::AssociationMap<ticl::mapWithSharedEnergy> layerClusterToAssociatedTracksterMap(layerClustersIds.size());
0258         std::vector<unsigned int> associatedRecoTracksterIndices;
0259         for (unsigned int i = 0; i < layerClustersIds.size(); ++i) {
0260           unsigned int layerClusterId = layerClustersIds[i];
0261           const auto& layerCluster = layerClusters[layerClusterId];
0262           float simFraction = 1.f / simTrackster.vertex_multiplicity(i);
0263           float squaredSimFraction = simFraction * simFraction;
0264           float squaredLayerClusterEnergy = layerCluster.energy() * layerCluster.energy();
0265           simToRecoScoresDenominator += squaredLayerClusterEnergy * squaredSimFraction;
0266           const auto& recoTracksterVec = layerClusterToTracksterMap[layerClusterId];
0267           for (const auto& recoTracksterElement : recoTracksterVec) {
0268             auto recoTracksterIndex = recoTracksterElement.index();
0269             auto recoSharedEnergy = recoTracksterElement.sharedEnergy();
0270             layerClusterToAssociatedTracksterMap[i].emplace_back(recoTracksterIndex, recoSharedEnergy);
0271             associatedRecoTracksterIndices.push_back(recoTracksterIndex);
0272           }
0273         }
0274         // keep only unique associatedRecoTracksterIndices
0275         std::sort(associatedRecoTracksterIndices.begin(), associatedRecoTracksterIndices.end());
0276         associatedRecoTracksterIndices.erase(
0277             std::unique(associatedRecoTracksterIndices.begin(), associatedRecoTracksterIndices.end()),
0278             associatedRecoTracksterIndices.end());
0279         // for each layer cluster, loop over associatedRecoTracksterIndices and add the missing reco tracksters with 0 shared energy to layerClusterToAssociatedTracksterMap
0280         for (unsigned int i = 0; i < layerClustersIds.size(); ++i) {
0281           unsigned int layerClusterId = layerClustersIds[i];
0282           const auto& recoTracksterVec = layerClusterToTracksterMap[layerClusterId];
0283           for (unsigned int recoTracksterIndex : associatedRecoTracksterIndices) {
0284             if (std::find_if(recoTracksterVec.begin(), recoTracksterVec.end(), [recoTracksterIndex](const auto& pair) {
0285                   return pair.index() == recoTracksterIndex;
0286                 }) == recoTracksterVec.end()) {
0287               layerClusterToAssociatedTracksterMap[i].emplace_back(recoTracksterIndex, 0.f);
0288             }
0289           }
0290         }
0291         assert(simToRecoScoresDenominator > 0.f);
0292         const float invDenominator = 1.f / simToRecoScoresDenominator;
0293 
0294         for (unsigned int i = 0; i < layerClustersIds.size(); ++i) {
0295           unsigned int layerClusterId = layerClustersIds[i];
0296           const auto& layerCluster = layerClusters[layerClusterId];
0297           float simFraction = 1.f / simTrackster.vertex_multiplicity(i);
0298           float squaredLayerClusterEnergy = layerCluster.energy() * layerCluster.energy();
0299           float simSharedEnergy = layerCluster.energy() * simFraction;
0300           float invLayerClusterEnergy = 1.f / layerCluster.energy();
0301           const auto& recoTracksterVec = layerClusterToAssociatedTracksterMap[i];
0302           for (const auto& recoTracksterElement : recoTracksterVec) {
0303             auto recoTracksterIndex = recoTracksterElement.index();
0304             float recoSharedEnergy = recoTracksterElement.sharedEnergy();
0305             edm::Ref<std::vector<ticl::Trackster>> recoTracksterRef(recoTrackstersHandle, recoTracksterIndex);
0306             float sharedEnergy = std::min(recoSharedEnergy, simSharedEnergy);
0307             float recoFraction =
0308                 recoSharedEnergy *
0309                 invLayerClusterEnergy;  // Either zero or one when no sharing of rechits between tracksters
0310             /* SimToReco score logic:
0311              - simFraction = 0 && recoFraction >= 0 : trackster contains non-sim associated elements : ignore in simToReco
0312              - simFraction > 0 && recoFraction == 0 : simhits not present in reco trackster : penalty in score by simFraction*E
0313              - simFraction == 1 && recoFraction == 1 : good association
0314              - 1 > simFraction > recoFraction > 0 :  we are missing some sim energy, penalty in score by the missing part : (simFraction-recoFraction)*E
0315              - 1 > recoFraction > simFraction > 0 : consider as good association, as there is enough reco to cover the sim
0316             */
0317             float simMinusRecoFraction = std::max(0.f, simFraction - recoFraction);
0318             float score = invDenominator * simMinusRecoFraction * simMinusRecoFraction * squaredLayerClusterEnergy;
0319             simTracksterToTracksterMap->insert(tracksterIndex, recoTracksterIndex, sharedEnergy, score);
0320           }
0321         }
0322       }
0323       auto sortingFunc = [](const auto& a, const auto& b) {
0324         if (a.score() != b.score())
0325           return a.score() < b.score();
0326         else
0327           return a.index() < b.index();
0328       };
0329 
0330       // Sort the maps by score in ascending order
0331       tracksterToSimTracksterMap->sort(sortingFunc);
0332       simTracksterToTracksterMap->sort(sortingFunc);
0333 
0334       // After populating the maps, store them in the event
0335       iEvent.put(std::move(tracksterToSimTracksterMap), tracksterToken.first + "To" + simTracksterToken.first);
0336       iEvent.put(std::move(simTracksterToTracksterMap), simTracksterToken.first + "To" + tracksterToken.first);
0337     }
0338   }
0339 }
0340 
0341 void AllTracksterToSimTracksterAssociatorsByLCsProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0342   edm::ParameterSetDescription desc;
0343   desc.add<std::string>("allLCtoTSAccoc", "allLayerClusterToTracksterAssociations");
0344   desc.add<std::vector<edm::InputTag>>(
0345       "tracksterCollections", {edm::InputTag("ticlTrackstersCLUE3DHigh"), edm::InputTag("ticlTrackstersLinks")});
0346   desc.add<std::vector<edm::InputTag>>(
0347       "simTracksterCollections", {edm::InputTag("ticlSimTracksters"), edm::InputTag("ticlSimTracksters", "fromCPs")});
0348   desc.add<edm::InputTag>("layerClusters", edm::InputTag("hgcalMergeLayerClusters"));
0349   descriptions.add("AllTracksterToSimTracksterAssociatorsByLCsProducer", desc);
0350 }
0351 
0352 // Define this as a plug-in
0353 DEFINE_FWK_MODULE(AllTracksterToSimTracksterAssociatorsByLCsProducer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!