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 "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
0013 #include "CommonTools/RecoAlgos/interface/MultiVectorManager.h"
0014 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
0015 #include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
0016 #include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
0017 
0018 class AllTracksterToSimTracksterAssociatorsByHitsProducer : public edm::global::EDProducer<> {
0019 public:
0020   explicit AllTracksterToSimTracksterAssociatorsByHitsProducer(const edm::ParameterSet&);
0021   ~AllTracksterToSimTracksterAssociatorsByHitsProducer() override = default;
0022 
0023   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0024 
0025 private:
0026   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0027 
0028   std::vector<std::pair<std::string, edm::EDGetTokenT<std::vector<ticl::Trackster>>>> tracksterCollectionTokens_;
0029   std::vector<std::pair<std::string, edm::EDGetTokenT<std::vector<ticl::Trackster>>>> simTracksterCollectionTokens_;
0030   std::vector<std::pair<std::string, edm::EDGetTokenT<ticl::AssociationMap<ticl::mapWithFraction>>>>
0031       hitToTracksterMapTokens_;
0032   std::vector<std::pair<std::string, edm::EDGetTokenT<ticl::AssociationMap<ticl::mapWithFraction>>>>
0033       tracksterToHitMapTokens_;
0034 
0035   std::vector<std::pair<std::string, edm::EDGetTokenT<ticl::AssociationMap<ticl::mapWithFraction>>>>
0036       hitToSimTracksterMapTokens_;
0037   std::vector<std::pair<std::string, edm::EDGetTokenT<ticl::AssociationMap<ticl::mapWithFraction>>>>
0038       simTracksterToHitMapTokens_;
0039 
0040   std::vector<edm::EDGetTokenT<HGCRecHitCollection>> hitsTokens_;
0041   edm::EDGetTokenT<std::vector<CaloParticle>> caloParticleToken_;
0042   edm::EDGetTokenT<ticl::AssociationMap<ticl::mapWithFraction>> hitToSimClusterMapToken_;
0043   edm::EDGetTokenT<ticl::AssociationMap<ticl::mapWithFraction>> hitToCaloParticleMapToken_;
0044 };
0045 
0046 AllTracksterToSimTracksterAssociatorsByHitsProducer::AllTracksterToSimTracksterAssociatorsByHitsProducer(
0047     const edm::ParameterSet& pset)
0048     : caloParticleToken_(consumes<std::vector<CaloParticle>>(pset.getParameter<edm::InputTag>("caloParticles"))),
0049       hitToSimClusterMapToken_(consumes<ticl::AssociationMap<ticl::mapWithFraction>>(
0050           pset.getParameter<edm::InputTag>("hitToSimClusterMap"))),
0051       hitToCaloParticleMapToken_(consumes<ticl::AssociationMap<ticl::mapWithFraction>>(
0052           pset.getParameter<edm::InputTag>("hitToCaloParticleMap"))) {
0053   const auto& tracksterCollections = pset.getParameter<std::vector<edm::InputTag>>("tracksterCollections");
0054 
0055   std::string allHitToTSAccoc = pset.getParameter<std::string>("allHitToTSAccoc");
0056   for (const auto& tag : tracksterCollections) {
0057     std::string label = tag.label();
0058     if (!tag.instance().empty()) {
0059       label += tag.instance();
0060     }
0061     tracksterCollectionTokens_.emplace_back(label, consumes<std::vector<ticl::Trackster>>(tag));
0062     hitToTracksterMapTokens_.emplace_back(
0063         label, consumes<ticl::AssociationMap<ticl::mapWithFraction>>(edm::InputTag(allHitToTSAccoc, "hitTo" + label)));
0064     tracksterToHitMapTokens_.emplace_back(
0065         label, consumes<ticl::AssociationMap<ticl::mapWithFraction>>(edm::InputTag(allHitToTSAccoc, label + "ToHit")));
0066   }
0067 
0068   const auto& simTracksterCollections = pset.getParameter<std::vector<edm::InputTag>>("simTracksterCollections");
0069   for (const auto& tag : simTracksterCollections) {
0070     std::string label = tag.label();
0071     if (!tag.instance().empty()) {
0072       label += tag.instance();
0073     }
0074     simTracksterCollectionTokens_.emplace_back(label, consumes<std::vector<ticl::Trackster>>(tag));
0075     hitToSimTracksterMapTokens_.emplace_back(
0076         label, consumes<ticl::AssociationMap<ticl::mapWithFraction>>(edm::InputTag(allHitToTSAccoc, "hitTo" + label)));
0077     simTracksterToHitMapTokens_.emplace_back(
0078         label, consumes<ticl::AssociationMap<ticl::mapWithFraction>>(edm::InputTag(allHitToTSAccoc, label + "ToHit")));
0079   }
0080 
0081   // Hits
0082   auto hitsTags = pset.getParameter<std::vector<edm::InputTag>>("hits");
0083   for (const auto& tag : hitsTags) {
0084     hitsTokens_.push_back(consumes<HGCRecHitCollection>(tag));
0085   }
0086 
0087   // Produce separate association maps for each trackster-simTrackster combination
0088   for (const auto& tracksterToken : tracksterCollectionTokens_) {
0089     for (const auto& simTracksterToken : simTracksterCollectionTokens_) {
0090       std::string instanceLabel = tracksterToken.first + "To" + simTracksterToken.first;
0091       produces<ticl::AssociationMap<ticl::mapWithSharedEnergyAndScore,
0092                                     std::vector<ticl::Trackster>,
0093                                     std::vector<ticl::Trackster>>>(instanceLabel);
0094       std::string reverseInstanceLabel = simTracksterToken.first + "To" + tracksterToken.first;
0095       produces<ticl::AssociationMap<ticl::mapWithSharedEnergyAndScore,
0096                                     std::vector<ticl::Trackster>,
0097                                     std::vector<ticl::Trackster>>>(reverseInstanceLabel);
0098     }
0099   }
0100 }
0101 
0102 void AllTracksterToSimTracksterAssociatorsByHitsProducer::produce(edm::StreamID,
0103                                                                   edm::Event& iEvent,
0104                                                                   const edm::EventSetup&) const {
0105   using namespace edm;
0106 
0107   MultiVectorManager<HGCRecHit> rechitManager;
0108   for (const auto& token : hitsTokens_) {
0109     Handle<HGCRecHitCollection> hitsHandle;
0110     iEvent.getByToken(token, hitsHandle);
0111 
0112     if (!hitsHandle.isValid()) {
0113       edm::LogWarning("AllTracksterToSimTracksterAssociatorsByHitsProducer")
0114           << "Missing HGCRecHitCollection for one of the hitsTokens.";
0115       continue;
0116     }
0117     rechitManager.addVector(*hitsHandle);
0118   }
0119 
0120   // Check if rechitManager is empty
0121   if (rechitManager.size() == 0) {
0122     edm::LogWarning("AllTracksterToSimTracksterAssociatorsByHitsProducer")
0123         << "No valid HGCRecHitCollections found. Association maps will be empty.";
0124 
0125     for (const auto& tracksterToken : tracksterCollectionTokens_) {
0126       for (const auto& simTracksterToken : simTracksterCollectionTokens_) {
0127         iEvent.put(std::make_unique<ticl::AssociationMap<ticl::mapWithSharedEnergyAndScore,
0128                                                          std::vector<ticl::Trackster>,
0129                                                          std::vector<ticl::Trackster>>>(),
0130                    tracksterToken.first + "To" + simTracksterToken.first);
0131         iEvent.put(std::make_unique<ticl::AssociationMap<ticl::mapWithSharedEnergyAndScore,
0132                                                          std::vector<ticl::Trackster>,
0133                                                          std::vector<ticl::Trackster>>>(),
0134                    simTracksterToken.first + "To" + tracksterToken.first);
0135       }
0136     }
0137     return;
0138   }
0139 
0140   Handle<ticl::AssociationMap<ticl::mapWithFraction>> hitToSimClusterMapHandle;
0141   iEvent.getByToken(hitToSimClusterMapToken_, hitToSimClusterMapHandle);
0142   const auto& hitToSimClusterMap = *hitToSimClusterMapHandle;
0143 
0144   Handle<ticl::AssociationMap<ticl::mapWithFraction>> hitToCaloParticleMapHandle;
0145   iEvent.getByToken(hitToCaloParticleMapToken_, hitToCaloParticleMapHandle);
0146   const auto& hitToCaloParticleMap = *hitToCaloParticleMapHandle;
0147 
0148   Handle<std::vector<CaloParticle>> caloParticlesHandle;
0149   iEvent.getByToken(caloParticleToken_, caloParticlesHandle);
0150 
0151   for (const auto& tracksterToken : tracksterCollectionTokens_) {
0152     Handle<std::vector<ticl::Trackster>> recoTrackstersHandle;
0153     iEvent.getByToken(tracksterToken.second, recoTrackstersHandle);
0154 
0155     if (!recoTrackstersHandle.isValid()) {
0156       edm::LogWarning("AllTracksterToSimTracksterAssociatorsByHitsProducer")
0157           << "No valid Trackster collection found. Association maps will be empty.";
0158       for (const auto& simTracksterToken : simTracksterCollectionTokens_) {
0159         Handle<std::vector<ticl::Trackster>> simTrackstersHandle;
0160 
0161         iEvent.put(std::make_unique<ticl::AssociationMap<ticl::mapWithSharedEnergyAndScore,
0162                                                          std::vector<ticl::Trackster>,
0163                                                          std::vector<ticl::Trackster>>>(),
0164                    tracksterToken.first + "To" + simTracksterToken.first);
0165         iEvent.put(std::make_unique<ticl::AssociationMap<ticl::mapWithSharedEnergyAndScore,
0166                                                          std::vector<ticl::Trackster>,
0167                                                          std::vector<ticl::Trackster>>>(),
0168                    simTracksterToken.first + "To" + tracksterToken.first);
0169       }
0170       continue;
0171     }
0172 
0173     const auto& recoTracksters = *recoTrackstersHandle;
0174 
0175     // Retrieve the correct HitToTracksterMap for the current trackster collection
0176     Handle<ticl::AssociationMap<ticl::mapWithFraction>> hitToTracksterMapHandle;
0177     auto tracksterMapTokenIter = std::find_if(
0178         hitToTracksterMapTokens_.begin(), hitToTracksterMapTokens_.end(), [&tracksterToken](const auto& pair) {
0179           return pair.first == tracksterToken.first;
0180         });
0181     if (tracksterMapTokenIter != hitToTracksterMapTokens_.end()) {
0182       iEvent.getByToken(tracksterMapTokenIter->second, hitToTracksterMapHandle);
0183     }
0184     const auto& hitToTracksterMap = *hitToTracksterMapHandle;
0185 
0186     // Retrieve the correct TracksterToHitMap for the current trackster collection
0187     Handle<ticl::AssociationMap<ticl::mapWithFraction>> tracksterToHitMapHandle;
0188     auto tracksterToHitMapTokenIter = std::find_if(
0189         tracksterToHitMapTokens_.begin(), tracksterToHitMapTokens_.end(), [&tracksterToken](const auto& pair) {
0190           return pair.first == tracksterToken.first;
0191         });
0192     if (tracksterToHitMapTokenIter != tracksterToHitMapTokens_.end()) {
0193       iEvent.getByToken(tracksterToHitMapTokenIter->second, tracksterToHitMapHandle);
0194     }
0195 
0196     if (!tracksterToHitMapHandle.isValid()) {
0197       edm::LogError("AllTracksterToSimTracksterAssociatorsByHitsProducer") << "tracksterToHitMapHandle is invalid";
0198       continue;
0199     }
0200 
0201     const auto& tracksterToHitMap = *tracksterToHitMapHandle;
0202 
0203     for (const auto& simTracksterToken : simTracksterCollectionTokens_) {
0204       Handle<std::vector<ticl::Trackster>> simTrackstersHandle;
0205       iEvent.getByToken(simTracksterToken.second, simTrackstersHandle);
0206 
0207       if (!simTrackstersHandle.isValid()) {
0208         iEvent.put(std::make_unique<ticl::AssociationMap<ticl::mapWithSharedEnergyAndScore,
0209                                                          std::vector<ticl::Trackster>,
0210                                                          std::vector<ticl::Trackster>>>(),
0211                    tracksterToken.first + "To" + simTracksterToken.first);
0212         iEvent.put(std::make_unique<ticl::AssociationMap<ticl::mapWithSharedEnergyAndScore,
0213                                                          std::vector<ticl::Trackster>,
0214                                                          std::vector<ticl::Trackster>>>(),
0215                    simTracksterToken.first + "To" + tracksterToken.first);
0216         continue;
0217       }
0218 
0219       const auto& simTracksters = *simTrackstersHandle;
0220 
0221       // Retrieve the correct HitToSimTracksterMap for the current simTrackster collection
0222       Handle<ticl::AssociationMap<ticl::mapWithFraction>> hitToSimTracksterMapHandle;
0223       auto simTracksterMapTokenIter =
0224           std::find_if(hitToSimTracksterMapTokens_.begin(),
0225                        hitToSimTracksterMapTokens_.end(),
0226                        [&simTracksterToken](const auto& pair) { return pair.first == simTracksterToken.first; });
0227       if (simTracksterMapTokenIter != hitToSimTracksterMapTokens_.end()) {
0228         iEvent.getByToken(simTracksterMapTokenIter->second, hitToSimTracksterMapHandle);
0229       }
0230       const auto& hitToSimTracksterMap = *hitToSimTracksterMapHandle;
0231 
0232       // Retrieve the correct SimTracksterToHitMap for the current simTrackster collection
0233       Handle<ticl::AssociationMap<ticl::mapWithFraction>> simTracksterToHitMapHandle;
0234       auto simTracksterToHitMapTokenIter =
0235           std::find_if(simTracksterToHitMapTokens_.begin(),
0236                        simTracksterToHitMapTokens_.end(),
0237                        [&simTracksterToken](const auto& pair) { return pair.first == simTracksterToken.first; });
0238       if (simTracksterToHitMapTokenIter != simTracksterToHitMapTokens_.end()) {
0239         iEvent.getByToken(simTracksterToHitMapTokenIter->second, simTracksterToHitMapHandle);
0240       }
0241       const auto& simTracksterToHitMap = *simTracksterToHitMapHandle;
0242 
0243       // Create the association maps
0244       auto tracksterToSimTracksterMap = std::make_unique<ticl::AssociationMap<ticl::mapWithSharedEnergyAndScore,
0245                                                                               std::vector<ticl::Trackster>,
0246                                                                               std::vector<ticl::Trackster>>>(
0247           recoTrackstersHandle, simTrackstersHandle, iEvent);
0248       auto simTracksterToTracksterMap = std::make_unique<ticl::AssociationMap<ticl::mapWithSharedEnergyAndScore,
0249                                                                               std::vector<ticl::Trackster>,
0250                                                                               std::vector<ticl::Trackster>>>(
0251           simTrackstersHandle, recoTrackstersHandle, iEvent);
0252 
0253       for (unsigned int tracksterIndex = 0; tracksterIndex < recoTracksters.size(); ++tracksterIndex) {
0254         edm::Ref<std::vector<ticl::Trackster>> recoTracksterRef(recoTrackstersHandle, tracksterIndex);
0255 
0256         float recoToSimScoresDenominator = 0.f;
0257         const auto& recoTracksterHitsAndFractions = tracksterToHitMap[tracksterIndex];
0258 
0259         if (tracksterToHitMap.size() == 0)
0260           continue;
0261 
0262         ticl::AssociationMap<ticl::mapWithFraction> hitToAssociatedSimTracksterMap(
0263             recoTracksterHitsAndFractions.size());
0264         std::vector<unsigned int> associatedSimTracksterIndices;
0265 
0266         for (unsigned int i = 0; i < recoTracksterHitsAndFractions.size(); ++i) {
0267           const auto& hitElement = recoTracksterHitsAndFractions[i];
0268           unsigned int hitIndex = hitElement.index();
0269           float recoFraction = hitElement.fraction();
0270           const auto& recHit = rechitManager[hitIndex];
0271           float squaredRecoFraction = recoFraction * recoFraction;
0272           float rechitEnergy = recHit.energy();
0273           float squaredRecHitEnergy = rechitEnergy * rechitEnergy;
0274           recoToSimScoresDenominator += squaredRecoFraction * squaredRecHitEnergy;
0275 
0276           const auto& hitToSimTracksterVec = hitToSimTracksterMap[hitIndex];
0277           for (const auto& simTracksterElement : hitToSimTracksterVec) {
0278             auto simTracksterIndex = simTracksterElement.index();
0279             const auto& simTrackster = simTracksters[simTracksterIndex];
0280             auto& seed = simTrackster.seedID();
0281             float simFraction = 0;
0282 
0283             if (seed == caloParticlesHandle.id()) {
0284               unsigned int caloParticleIndex = simTrackster.seedIndex();
0285               auto it =
0286                   std::find_if(hitToCaloParticleMap[hitIndex].begin(),
0287                                hitToCaloParticleMap[hitIndex].end(),
0288                                [caloParticleIndex](const auto& pair) { return pair.index() == caloParticleIndex; });
0289               if (it != hitToCaloParticleMap[hitIndex].end()) {
0290                 simFraction = it->fraction();
0291               }
0292             } else {
0293               unsigned int simClusterIndex = simTracksters[simTracksterIndex].seedIndex();
0294               auto it = std::find_if(hitToSimClusterMap[hitIndex].begin(),
0295                                      hitToSimClusterMap[hitIndex].end(),
0296                                      [simClusterIndex](const auto& pair) { return pair.index() == simClusterIndex; });
0297               if (it != hitToSimClusterMap[hitIndex].end()) {
0298                 simFraction = it->fraction();
0299               }
0300             }
0301 
0302             hitToAssociatedSimTracksterMap.insert(i, simTracksterIndex, simFraction);
0303             associatedSimTracksterIndices.push_back(simTracksterIndex);
0304           }
0305         }
0306         std::sort(associatedSimTracksterIndices.begin(), associatedSimTracksterIndices.end());
0307         associatedSimTracksterIndices.erase(
0308             std::unique(associatedSimTracksterIndices.begin(), associatedSimTracksterIndices.end()),
0309             associatedSimTracksterIndices.end());
0310 
0311         // Add missing sim tracksters with 0 shared energy to hitToAssociatedSimTracksterMap
0312         for (unsigned int i = 0; i < recoTracksterHitsAndFractions.size(); ++i) {
0313           unsigned int hitId = recoTracksterHitsAndFractions[i].index();
0314           const auto& simTracksterVec = hitToSimTracksterMap[hitId];
0315           for (unsigned int simTracksterIndex : associatedSimTracksterIndices) {
0316             if (std::find_if(simTracksterVec.begin(), simTracksterVec.end(), [simTracksterIndex](const auto& pair) {
0317                   return pair.index() == simTracksterIndex;
0318                 }) == simTracksterVec.end()) {
0319               hitToAssociatedSimTracksterMap.insert(i, simTracksterIndex, 0);
0320             }
0321           }
0322         }
0323 
0324         const float invDenominator = 1.f / recoToSimScoresDenominator;
0325 
0326         for (unsigned int i = 0; i < recoTracksterHitsAndFractions.size(); ++i) {
0327           unsigned int hitIndex = recoTracksterHitsAndFractions[i].index();
0328           const auto& recHit = rechitManager[hitIndex];
0329           float recoFraction = recoTracksterHitsAndFractions[i].fraction();
0330           float rechitEnergy = recHit.energy();
0331           float squaredRecHitEnergy = rechitEnergy * rechitEnergy;
0332           float recoSharedEnergy = rechitEnergy * recoFraction;
0333           const auto& simTracksterVec = hitToAssociatedSimTracksterMap[i];
0334           for (const auto& simTracksterElement : simTracksterVec) {
0335             auto simTracksterIndex = simTracksterElement.index();
0336             auto simFraction = simTracksterElement.fraction();
0337             edm::Ref<std::vector<ticl::Trackster>> simTracksterRef(simTrackstersHandle, simTracksterIndex);
0338             float sharedEnergy = std::min(simFraction * rechitEnergy, recoSharedEnergy);
0339             /* RecoToSim score logic:
0340              - simFraction >= 0 && recoFraction == 0 : simtrackster contains non-reco associated elements : ignore in recoToSim association
0341              - simFraction == 0 && recoFraction > 0 : rechits not present in sim trackster : penalty in score by recoFraction*E
0342              - simFraction == 1 && recoFraction == 1 : good association
0343              - 1 > recoFraction > simFraction > 0 :  sim does not contain some reco energy, penalty in score by the additional part : (recoFraction-simFraction)*E
0344              - 1 > simFraction> recoFraction > 0 : consider as good association, as there is enough sim to cover the reco
0345             */
0346             float recoMinusSimFraction = std::max(0.f, recoFraction - simFraction);
0347             float score = invDenominator * recoMinusSimFraction * recoMinusSimFraction * squaredRecHitEnergy;
0348             tracksterToSimTracksterMap->insert(recoTracksterRef, simTracksterRef, sharedEnergy, score);
0349           }
0350         }
0351       }
0352 
0353       // Reverse mapping: SimTrackster -> RecoTrackster
0354       for (unsigned int tracksterIndex = 0; tracksterIndex < simTracksters.size(); ++tracksterIndex) {
0355         edm::Ref<std::vector<ticl::Trackster>> simTracksterRef(simTrackstersHandle, tracksterIndex);
0356         float simToRecoScoresDenominator = 0.f;
0357         const auto& simTracksterHitsAndFractions = simTracksterToHitMap[tracksterIndex];
0358         ticl::AssociationMap<ticl::mapWithFraction> hitToAssociatedRecoTracksterMap(
0359             simTracksterHitsAndFractions.size());
0360         std::vector<unsigned int> associatedRecoTracksterIndices;
0361         const auto& simTrackster = simTracksters[tracksterIndex];
0362         auto& seed = simTrackster.seedID();
0363         unsigned int simObjectIndex = simTrackster.seedIndex();
0364         bool isSimTracksterFromCP = (seed == caloParticlesHandle.id());
0365         std::vector<float> simFractions(simTracksterHitsAndFractions.size(), 0.f);
0366         for (unsigned int i = 0; i < simTracksterHitsAndFractions.size(); ++i) {
0367           auto hitIndex = simTracksterHitsAndFractions[i].index();
0368           auto it = isSimTracksterFromCP
0369                         ? (std::find_if(hitToCaloParticleMap[hitIndex].begin(),
0370                                         hitToCaloParticleMap[hitIndex].end(),
0371                                         [simObjectIndex](const auto& pair) { return pair.index() == simObjectIndex; }))
0372                         : std::find_if(hitToSimClusterMap[hitIndex].begin(),
0373                                        hitToSimClusterMap[hitIndex].end(),
0374                                        [simObjectIndex](const auto& pair) { return pair.index() == simObjectIndex; });
0375           if ((isSimTracksterFromCP and it != hitToCaloParticleMap[hitIndex].end()) or
0376               (!isSimTracksterFromCP and it != hitToSimClusterMap[hitIndex].end())) {
0377             simFractions[i] = it->fraction();
0378           }
0379           float simFraction = simFractions[i];
0380           const auto& recHit = rechitManager[hitIndex];
0381           float rechitEnergy = recHit.energy();
0382           float squaredSimFraction = simFraction * simFraction;
0383           float squaredRecHitEnergy = rechitEnergy * rechitEnergy;
0384           simToRecoScoresDenominator += squaredSimFraction * squaredRecHitEnergy;
0385           const auto& hitToRecoTracksterVec = hitToTracksterMap[hitIndex];
0386           for (const auto& recoTracksterElement : hitToRecoTracksterVec) {
0387             unsigned int recoTracksterIndex = recoTracksterElement.index();
0388             float recoFraction = recoTracksterElement.fraction();
0389             hitToAssociatedRecoTracksterMap.insert(i, recoTracksterIndex, recoFraction);
0390             associatedRecoTracksterIndices.push_back(recoTracksterIndex);
0391           }
0392         }
0393 
0394         std::sort(associatedRecoTracksterIndices.begin(), associatedRecoTracksterIndices.end());
0395         associatedRecoTracksterIndices.erase(
0396             std::unique(associatedRecoTracksterIndices.begin(), associatedRecoTracksterIndices.end()),
0397             associatedRecoTracksterIndices.end());
0398 
0399         for (unsigned int i = 0; i < simTracksterHitsAndFractions.size(); ++i) {
0400           unsigned int hitIndex = simTracksterHitsAndFractions[i].index();
0401           const auto& hitToRecoTracksterVec = hitToTracksterMap[hitIndex];
0402           for (unsigned int recoTracksterIndex : associatedRecoTracksterIndices) {
0403             if (std::find_if(
0404                     hitToRecoTracksterVec.begin(), hitToRecoTracksterVec.end(), [recoTracksterIndex](const auto& pair) {
0405                       return pair.index() == recoTracksterIndex;
0406                     }) == hitToRecoTracksterVec.end()) {
0407               hitToAssociatedRecoTracksterMap.insert(i, recoTracksterIndex, 0);
0408             }
0409           }
0410         }
0411 
0412         assert(simToRecoScoresDenominator > 0.f);
0413         const float invDenominator = 1.f / simToRecoScoresDenominator;
0414         for (unsigned int i = 0; i < simTracksterHitsAndFractions.size(); ++i) {
0415           const auto& hitIndex = simTracksterHitsAndFractions[i].index();
0416           float simFraction = simFractions[i];
0417           const auto& recHit = rechitManager[hitIndex];
0418           float rechitEnergy = recHit.energy();
0419           float squaredRecHitEnergy = rechitEnergy * rechitEnergy;
0420           float simSharedEnergy = rechitEnergy * simFraction;
0421           const auto& hitToRecoTracksterVec = hitToAssociatedRecoTracksterMap[i];
0422           for (const auto& recoTracksterElement : hitToRecoTracksterVec) {
0423             auto recoTracksterIndex = recoTracksterElement.index();
0424             float recoFraction =
0425                 recoTracksterElement.fraction();  // Either zero or one when no sharing of rechits between tracksters
0426             edm::Ref<std::vector<ticl::Trackster>> recoTracksterRef(recoTrackstersHandle, recoTracksterIndex);
0427             float sharedEnergy = std::min(recoFraction * rechitEnergy, simSharedEnergy);
0428             /* SimToReco score logic:
0429              - simFraction = 0 && recoFraction >= 0 : trackster contains non-sim associated elements : ignore in simToReco
0430              - simFraction > 0 && recoFraction == 0 : simhits not present in reco trackster : penalty in score by simFraction*E
0431              - simFraction == 1 && recoFraction == 1 : good association
0432              - 1 > simFraction > recoFraction > 0 :  we are missing some sim energy, penalty in score by the missing part : (simFraction-recoFraction)*E
0433              - 1 > recoFraction > simFraction > 0 : consider as good association, as there is enough reco to cover the sim
0434             */
0435             float simMinusRecoFraction = std::max(0.f, simFraction - recoFraction);
0436             float score = invDenominator * simMinusRecoFraction * simMinusRecoFraction * squaredRecHitEnergy;
0437             simTracksterToTracksterMap->insert(simTracksterRef, recoTracksterRef, sharedEnergy, score);
0438           }
0439         }
0440       }
0441 
0442       auto sortingFunc = [](const auto& a, const auto& b) {
0443         if (a.score() != b.score())
0444           return a.score() < b.score();
0445         else
0446           return a.index() < b.index();
0447       };
0448 
0449       // Sort the maps by score in ascending order
0450       tracksterToSimTracksterMap->sort(sortingFunc);
0451       simTracksterToTracksterMap->sort(sortingFunc);
0452 
0453       // After populating the maps, store them in the event
0454       iEvent.put(std::move(tracksterToSimTracksterMap), tracksterToken.first + "To" + simTracksterToken.first);
0455       iEvent.put(std::move(simTracksterToTracksterMap), simTracksterToken.first + "To" + tracksterToken.first);
0456     }
0457   }
0458 }
0459 
0460 void AllTracksterToSimTracksterAssociatorsByHitsProducer::fillDescriptions(
0461     edm::ConfigurationDescriptions& descriptions) {
0462   edm::ParameterSetDescription desc;
0463   desc.add<std::string>("allHitToTSAccoc", "allHitToTracksterAssociations");
0464   desc.add<std::vector<edm::InputTag>>(
0465       "tracksterCollections", {edm::InputTag("ticlTrackstersCLUE3DHigh"), edm::InputTag("ticlTrackstersLinks")});
0466   desc.add<std::vector<edm::InputTag>>(
0467       "simTracksterCollections", {edm::InputTag("ticlSimTracksters"), edm::InputTag("ticlSimTracksters", "fromCPs")});
0468   desc.add<std::vector<edm::InputTag>>("hits",
0469                                        {edm::InputTag("HGCalRecHit", "HGCEERecHits"),
0470                                         edm::InputTag("HGCalRecHit", "HGCHEFRecHits"),
0471                                         edm::InputTag("HGCalRecHit", "HGCHEBRecHits")});
0472   desc.add<edm::InputTag>("hitToSimClusterMap",
0473                           edm::InputTag("hitToSimClusterCaloParticleAssociator", "hitToSimClusterMap"));
0474   desc.add<edm::InputTag>("hitToCaloParticleMap",
0475                           edm::InputTag("hitToSimClusterCaloParticleAssociator", "hitToCaloParticleMap"));
0476   desc.add<edm::InputTag>("caloParticles", edm::InputTag("mix", "MergedCaloTruth"));
0477 
0478   descriptions.add("AllTracksterToSimTracksterAssociatorsByHitsProducer", desc);
0479 }
0480 
0481 // Define this as a plug-in
0482 DEFINE_FWK_MODULE(AllTracksterToSimTracksterAssociatorsByHitsProducer);
[ 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!