TICL/plugins/SimTrackstersProducer.cc

0001 // Author: Felice Pantaleo, Leonardo Cristella - felice.pantaleo@cern.ch, leonardo.cristella@cern.ch
0002 // Date: 09/2021
0003
0004 // user include files
0005
0006 #include "FWCore/Framework/interface/ESHandle.h"
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "FWCore/Framework/interface/Frameworkfwd.h"
0009 #include "FWCore/Framework/interface/MakerMacros.h"
0010 #include "FWCore/Framework/interface/stream/EDProducer.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0014 #include "DataFormats/Common/interface/OrphanHandle.h"
0015
0016 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
0017 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0018
0019 #include "DataFormats/HGCalReco/interface/Trackster.h"
0020 #include "DataFormats/HGCalReco/interface/TICLCandidate.h"
0021 #include "DataFormats/TrackReco/interface/Track.h"
0022 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0023 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0024
0025 #include "DataFormats/Common/interface/ValueMap.h"
0026 #include "SimDataFormats/Associations/interface/LayerClusterToSimClusterAssociator.h"
0027 #include "SimDataFormats/Associations/interface/LayerClusterToCaloParticleAssociator.h"
0028
0029 #include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
0030 #include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
0031 #include "SimDataFormats/CaloAnalysis/interface/MtdSimTrackster.h"
0032 #include "SimDataFormats/CaloAnalysis/interface/MtdSimTracksterFwd.h"
0033 #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
0034
0035 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0036 #include "SimDataFormats/TrackingAnalysis/interface/UniqueSimTrackId.h"
0037
0038 #include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h"
0039
0040 #include "DataFormats/HGCalReco/interface/Common.h"
0041
0042 #include <CLHEP/Units/SystemOfUnits.h>
0043
0044 #include "TrackstersPCA.h"
0045 #include <vector>
0046 #include <map>
0047 #include <iterator>
0048 #include <algorithm>
0049 #include <numeric>
0050
0051 using namespace ticl;
0052
0053 class SimTrackstersProducer : public edm::stream::EDProducer<> {
0054 public:
0055   explicit SimTrackstersProducer(const edm::ParameterSet&);
0056   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0057
0058   void produce(edm::Event&, const edm::EventSetup&) override;
0059   void makePUTrackster(const std::vector<float>& inputClusterMask,
0060                        std::vector<float>& output_mask,
0061                        std::vector<Trackster>& result,
0062                        const edm::ProductID seed,
0063                        int loop_index);
0064
0065   void addTrackster(const int index,
0066                     const std::vector<std::pair<edm::Ref<reco::CaloClusterCollection>, std::pair<float, float>>>& lcVec,
0067                     const std::vector<float>& inputClusterMask,
0068                     const float fractionCut_,
0069                     const float energy,
0070                     const int pdgId,
0071                     const int charge,
0072                     const float time,
0073                     const edm::ProductID seed,
0074                     const Trackster::IterationIndex iter,
0075                     std::vector<float>& output_mask,
0076                     std::vector<Trackster>& result,
0077                     int& loop_index,
0078                     const bool add = false);
0079
0080 private:
0081   std::string detector_;
0082   const bool doNose_ = false;
0083   const bool computeLocalTime_;
0084   const edm::EDGetTokenT<std::vector<reco::CaloCluster>> clusters_token_;
0085   const edm::EDGetTokenT<edm::ValueMap<std::pair<float, float>>> clustersTime_token_;
0086   const edm::EDGetTokenT<std::vector<float>> filtered_layerclusters_mask_token_;
0087
0088   const edm::EDGetTokenT<std::vector<SimCluster>> simclusters_token_;
0089   const edm::EDGetTokenT<std::vector<CaloParticle>> caloparticles_token_;
0090   const edm::EDGetTokenT<MtdSimTracksterCollection> MTDSimTrackstersToken_;
0091
0092   const edm::EDGetTokenT<ticl::SimToRecoCollectionWithSimClusters> associatorMapSimClusterToReco_token_;
0093   const edm::EDGetTokenT<ticl::SimToRecoCollection> associatorMapCaloParticleToReco_token_;
0094   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geom_token_;
0095   hgcal::RecHitTools rhtools_;
0096   const float fractionCut_;
0097   const float qualityCutTrack_;
0098   const edm::EDGetTokenT<std::vector<TrackingParticle>> trackingParticleToken_;
0099
0100   const edm::EDGetTokenT<std::vector<reco::Track>> recoTracksToken_;
0101   const StringCutObjectSelector<reco::Track> cutTk_;
0102
0103   const edm::EDGetTokenT<reco::SimToRecoCollection> associatormapStRsToken_;
0104   const edm::EDGetTokenT<reco::RecoToSimCollection> associatormapRtSsToken_;
0105   const edm::EDGetTokenT<SimTrackToTPMap> associationSimTrackToTPToken_;
0106 };
0107
0108 DEFINE_FWK_MODULE(SimTrackstersProducer);
0109
0110 SimTrackstersProducer::SimTrackstersProducer(const edm::ParameterSet& ps)
0111     : detector_(ps.getParameter<std::string>("detector")),
0112       doNose_(detector_ == "HFNose"),
0113       computeLocalTime_(ps.getParameter<bool>("computeLocalTime")),
0114       clusters_token_(consumes(ps.getParameter<edm::InputTag>("layer_clusters"))),
0115       clustersTime_token_(consumes(ps.getParameter<edm::InputTag>("time_layerclusters"))),
0116       filtered_layerclusters_mask_token_(consumes(ps.getParameter<edm::InputTag>("filtered_mask"))),
0117       simclusters_token_(consumes(ps.getParameter<edm::InputTag>("simclusters"))),
0118       caloparticles_token_(consumes(ps.getParameter<edm::InputTag>("caloparticles"))),
0119       MTDSimTrackstersToken_(consumes<MtdSimTracksterCollection>(ps.getParameter<edm::InputTag>("MtdSimTracksters"))),
0120       associatorMapSimClusterToReco_token_(
0121           consumes(ps.getParameter<edm::InputTag>("layerClusterSimClusterAssociator"))),
0122       associatorMapCaloParticleToReco_token_(
0123           consumes(ps.getParameter<edm::InputTag>("layerClusterCaloParticleAssociator"))),
0124       geom_token_(esConsumes()),
0125       fractionCut_(ps.getParameter<double>("fractionCut")),
0126       qualityCutTrack_(ps.getParameter<double>("qualityCutTrack")),
0127       trackingParticleToken_(
0128           consumes<std::vector<TrackingParticle>>(ps.getParameter<edm::InputTag>("trackingParticles"))),
0129       recoTracksToken_(consumes<std::vector<reco::Track>>(ps.getParameter<edm::InputTag>("recoTracks"))),
0130       cutTk_(ps.getParameter<std::string>("cutTk")),
0131       associatormapStRsToken_(consumes(ps.getParameter<edm::InputTag>("tpToTrack"))),
0132       associationSimTrackToTPToken_(consumes(ps.getParameter<edm::InputTag>("simTrackToTPMap"))) {
0133   produces<TracksterCollection>();
0134   produces<std::vector<float>>();
0135   produces<TracksterCollection>("fromCPs");
0136   produces<TracksterCollection>("PU");
0137   produces<std::vector<float>>("fromCPs");
0138   produces<std::map<uint, std::vector<uint>>>();
0139   produces<std::vector<TICLCandidate>>();
0140 }
0141
0142 void SimTrackstersProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0143   edm::ParameterSetDescription desc;
0144   desc.add<std::string>("detector", "HGCAL");
0145   desc.add<bool>("computeLocalTime", "false");
0146   desc.add<edm::InputTag>("layer_clusters", edm::InputTag("hgcalMergeLayerClusters"));
0147   desc.add<edm::InputTag>("time_layerclusters", edm::InputTag("hgcalMergeLayerClusters", "timeLayerCluster"));
0148   desc.add<edm::InputTag>("filtered_mask", edm::InputTag("filteredLayerClustersSimTracksters", "ticlSimTracksters"));
0149   desc.add<edm::InputTag>("simclusters", edm::InputTag("mix", "MergedCaloTruth"));
0150   desc.add<edm::InputTag>("caloparticles", edm::InputTag("mix", "MergedCaloTruth"));
0151   desc.add<edm::InputTag>("MtdSimTracksters", edm::InputTag("mix", "MergedMtdTruthST"));
0152   desc.add<edm::InputTag>("layerClusterSimClusterAssociator",
0153                           edm::InputTag("layerClusterSimClusterAssociationProducer"));
0154   desc.add<edm::InputTag>("layerClusterCaloParticleAssociator",
0155                           edm::InputTag("layerClusterCaloParticleAssociationProducer"));
0156   desc.add<edm::InputTag>("recoTracks", edm::InputTag("generalTracks"));
0157   desc.add<std::string>("cutTk",
0158                         "1.48 < abs(eta) < 3.0 && pt > 1. && quality(\"highPurity\") && "
0159                         "hitPattern().numberOfLostHits(\"MISSING_OUTER_HITS\") < 5");
0160   desc.add<edm::InputTag>("tpToTrack", edm::InputTag("trackingParticleRecoTrackAsssociation"));
0161
0162   desc.add<edm::InputTag>("trackingParticles", edm::InputTag("mix", "MergedTrackTruth"));
0163
0164   desc.add<edm::InputTag>("simTrackToTPMap", edm::InputTag("simHitTPAssocProducer", "simTrackToTP"));
0165   desc.add<double>("fractionCut", 0.);
0166   desc.add<double>("qualityCutTrack", 0.75);
0167
0168   descriptions.addWithDefaultLabel(desc);
0169 }
0170 void SimTrackstersProducer::makePUTrackster(const std::vector<float>& inputClusterMask,
0171                                             std::vector<float>& output_mask,
0172                                             std::vector<Trackster>& result,
0173                                             const edm::ProductID seed,
0174                                             int loop_index) {
0175   Trackster tmpTrackster;
0176   for (size_t i = 0; i < output_mask.size(); i++) {
0177     const float remaining_fraction = output_mask[i];
0178     if (remaining_fraction > std::numeric_limits<float>::epsilon()) {
0179       tmpTrackster.vertices().push_back(i);
0180       tmpTrackster.vertex_multiplicity().push_back(1. / remaining_fraction);
0181     }
0182   }
0183   tmpTrackster.setSeed(seed, 0);
0184   result.push_back(tmpTrackster);
0185 }
0186
0187 void SimTrackstersProducer::addTrackster(
0188     const int index,
0189     const std::vector<std::pair<edm::Ref<reco::CaloClusterCollection>, std::pair<float, float>>>& lcVec,
0190     const std::vector<float>& inputClusterMask,
0191     const float fractionCut_,
0192     const float energy,
0193     const int pdgId,
0194     const int charge,
0195     const float time,
0196     const edm::ProductID seed,
0197     const Trackster::IterationIndex iter,
0198     std::vector<float>& output_mask,
0199     std::vector<Trackster>& result,
0200     int& loop_index,
0201     const bool add) {
0202   Trackster tmpTrackster;
0203   if (lcVec.empty()) {
0204     result[index] = tmpTrackster;
0205     return;
0206   }
0207
0208   tmpTrackster.vertices().reserve(lcVec.size());
0209   tmpTrackster.vertex_multiplicity().reserve(lcVec.size());
0210   for (auto const& [lc, energyScorePair] : lcVec) {
0211     if (inputClusterMask[lc.index()] > 0) {
0212       float fraction = energyScorePair.first / lc->energy();
0213       if (fraction < fractionCut_)
0214         continue;
0215       tmpTrackster.vertices().push_back(lc.index());
0216       output_mask[lc.index()] -= fraction;
0217       tmpTrackster.vertex_multiplicity().push_back(1. / fraction);
0218     }
0219   }
0220
0221   tmpTrackster.setIdProbability(tracksterParticleTypeFromPdgId(pdgId, charge), 1.f);
0222   tmpTrackster.setRegressedEnergy(energy);
0223   tmpTrackster.setIteration(iter);
0224   tmpTrackster.setSeed(seed, index);
0225   tmpTrackster.setBoundaryTime(time);
0226   if (add) {
0227     result[index] = tmpTrackster;
0228     loop_index += 1;
0229   } else {
0230     result.push_back(tmpTrackster);
0231   }
0232 }
0233
0234 void SimTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0235   auto result = std::make_unique<TracksterCollection>();
0236   auto output_mask = std::make_unique<std::vector<float>>();
0237   auto result_fromCP = std::make_unique<TracksterCollection>();
0238   auto resultPU = std::make_unique<TracksterCollection>();
0239   auto output_mask_fromCP = std::make_unique<std::vector<float>>();
0240   auto cpToSc_SimTrackstersMap = std::make_unique<std::map<uint, std::vector<uint>>>();
0241   auto result_ticlCandidates = std::make_unique<std::vector<TICLCandidate>>();
0242
0243   const auto& layerClustersHandle = evt.getHandle(clusters_token_);
0244   const auto& layerClustersTimesHandle = evt.getHandle(clustersTime_token_);
0245   const auto& inputClusterMaskHandle = evt.getHandle(filtered_layerclusters_mask_token_);
0246
0247   // Validate input collections
0248   if (!layerClustersHandle.isValid() || !layerClustersTimesHandle.isValid() || !inputClusterMaskHandle.isValid()) {
0249     edm::LogWarning("SimTrackstersProducer") << "Missing input collections. Producing empty outputs.";
0250
0251     evt.put(std::move(result_ticlCandidates));
0252     evt.put(std::move(output_mask));
0253     evt.put(std::move(result_fromCP), "fromCPs");
0254     evt.put(std::move(resultPU), "PU");
0255     evt.put(std::move(output_mask_fromCP), "fromCPs");
0256     evt.put(std::move(cpToSc_SimTrackstersMap));
0257     return;
0258   }
0259
0260   // Proceed if inputs are valid
0261   const auto& layerClusters = *layerClustersHandle;
0262   const auto& layerClustersTimes = *layerClustersTimesHandle;
0263   const auto& inputClusterMask = *inputClusterMaskHandle;
0264
0265   output_mask->resize(layerClusters.size(), 1.f);
0266   output_mask_fromCP->resize(layerClusters.size(), 1.f);
0267
0268   const auto& simclusters = evt.get(simclusters_token_);
0269   edm::Handle<std::vector<CaloParticle>> caloParticles_h;
0270   evt.getByToken(caloparticles_token_, caloParticles_h);
0271   const auto& caloparticles = *caloParticles_h;
0272
0273   edm::Handle<MtdSimTracksterCollection> MTDSimTracksters_h;
0274   evt.getByToken(MTDSimTrackstersToken_, MTDSimTracksters_h);
0275
0276   const auto& simClustersToRecoColl = evt.get(associatorMapSimClusterToReco_token_);
0277   const auto& caloParticlesToRecoColl = evt.get(associatorMapCaloParticleToReco_token_);
0278
0279   edm::Handle<std::vector<TrackingParticle>> trackingParticles_h;
0280   evt.getByToken(trackingParticleToken_, trackingParticles_h);
0281   edm::Handle<std::vector<reco::Track>> recoTracks_h;
0282   evt.getByToken(recoTracksToken_, recoTracks_h);
0283   const auto& TPtoRecoTrackMap = evt.get(associatormapStRsToken_);
0284   const auto& simTrackToTPMap = evt.get(associationSimTrackToTPToken_);
0285   const auto& recoTracks = *recoTracks_h;
0286
0287   const auto& geom = es.getData(geom_token_);
0288   rhtools_.setGeometry(geom);
0289   const auto num_simclusters = simclusters.size();
0290   result->reserve(num_simclusters);  // Conservative size, will call shrink_to_fit later
0291   const auto num_caloparticles = caloparticles.size();
0292   result_fromCP->resize(num_caloparticles);
0293   std::map<uint, uint> SimClusterToCaloParticleMap;
0294   int loop_index = 0;
0295   for (const auto& [key, lcVec] : caloParticlesToRecoColl) {
0296     auto const& cp = *(key);
0297     auto cpIndex = &cp - &caloparticles[0];
0298     for (const auto& scRef : cp.simClusters()) {
0299       auto const& sc = *(scRef);
0300       auto const scIndex = &sc - &simclusters[0];
0301       SimClusterToCaloParticleMap[scIndex] = cpIndex;
0302     }
0303
0304     auto regr_energy = cp.energy();
0305     std::vector<uint> scSimTracksterIdx;
0306     scSimTracksterIdx.reserve(cp.simClusters().size());
0307
0308     // Create a Trackster from the object entering HGCal
0309     if (cp.g4Tracks()[0].crossedBoundary()) {
0310       regr_energy = cp.g4Tracks()[0].getMomentumAtBoundary().energy();
0311       float time = cp.g4Tracks()[0].getPositionAtBoundary().t() *
0312                    CLHEP::s;  // Geant4 time is in seconds, convert to ns (CLHEP::s = 1e9)
0313       addTrackster(cpIndex,
0314                    lcVec,
0315                    inputClusterMask,
0316                    fractionCut_,
0317                    regr_energy,
0318                    cp.pdgId(),
0319                    cp.charge(),
0320                    time,
0321                    key.id(),
0322                    ticl::Trackster::SIM,
0323                    *output_mask,
0324                    *result,
0325                    loop_index);
0326     } else {
0327       for (const auto& scRef : cp.simClusters()) {
0328         const auto& it = simClustersToRecoColl.find(scRef);
0329         if (it == simClustersToRecoColl.end())
0330           continue;
0331         const auto& lcVec = it->val;
0332         auto const& sc = *(scRef);
0333         auto const scIndex = &sc - &simclusters[0];
0334
0335         addTrackster(scIndex,
0336                      lcVec,
0337                      inputClusterMask,
0338                      fractionCut_,
0339                      sc.g4Tracks()[0].getMomentumAtBoundary().energy(),
0340                      sc.pdgId(),
0341                      sc.charge(),
0342                      sc.g4Tracks()[0].getPositionAtBoundary().t() *
0343                          CLHEP::s,  // Geant4 time is in seconds, convert to ns (CLHEP::s = 1e9)
0344                      scRef.id(),
0345                      ticl::Trackster::SIM,
0346                      *output_mask,
0347                      *result,
0348                      loop_index);
0349
0350         if (result->empty())
0351           continue;
0352         const auto index = result->size() - 1;
0353         if (std::find(scSimTracksterIdx.begin(), scSimTracksterIdx.end(), index) == scSimTracksterIdx.end()) {
0354           scSimTracksterIdx.emplace_back(index);
0355         }
0356       }
0357       scSimTracksterIdx.shrink_to_fit();
0358     }
0359     float time = cp.simTime();
0360     // Create a Trackster from any CP
0361     addTrackster(cpIndex,
0362                  lcVec,
0363                  inputClusterMask,
0364                  fractionCut_,
0365                  regr_energy,
0366                  cp.pdgId(),
0367                  cp.charge(),
0368                  time,
0369                  key.id(),
0370                  ticl::Trackster::SIM_CP,
0371                  *output_mask_fromCP,
0372                  *result_fromCP,
0373                  loop_index,
0374                  true);
0375
0376     if (result_fromCP->empty())
0377       continue;
0378     const auto index = loop_index - 1;
0379     if (cpToSc_SimTrackstersMap->find(index) == cpToSc_SimTrackstersMap->end()) {
0380       (*cpToSc_SimTrackstersMap)[index] = scSimTracksterIdx;
0381     }
0382   }
0383   // TODO: remove time computation from PCA calculation and
0384   //       store time from boundary position in simTracksters
0385   ticl::assignPCAtoTracksters(*result,
0386                               layerClusters,
0387                               layerClustersTimes,
0388                               rhtools_.getPositionLayer(rhtools_.lastLayerEE(doNose_)).z(),
0389                               rhtools_,
0390                               computeLocalTime_);
0391   result->shrink_to_fit();
0392   ticl::assignPCAtoTracksters(*result_fromCP,
0393                               layerClusters,
0394                               layerClustersTimes,
0395                               rhtools_.getPositionLayer(rhtools_.lastLayerEE(doNose_)).z(),
0396                               rhtools_,
0397                               computeLocalTime_);
0398
0399   makePUTrackster(inputClusterMask, *output_mask, *resultPU, caloParticles_h.id(), 0);
0400
0401   auto simTrackToRecoTrack = [&](UniqueSimTrackId simTkId) -> std::pair<int, float> {
0402     int trackIdx = -1;
0403     float quality = 0.f;
0404     auto ipos = simTrackToTPMap.mapping.find(simTkId);
0405     if (ipos != simTrackToTPMap.mapping.end()) {
0406       auto jpos = TPtoRecoTrackMap.find((ipos->second));
0407       if (jpos != TPtoRecoTrackMap.end()) {
0408         auto& associatedRecoTracks = jpos->val;
0409         if (!associatedRecoTracks.empty()) {
0410           // associated reco tracks are sorted by decreasing quality
0411           if (associatedRecoTracks[0].second > qualityCutTrack_) {
0412             trackIdx = &(*associatedRecoTracks[0].first) - &recoTracks[0];
0413             quality = associatedRecoTracks[0].second;
0414           }
0415         }
0416       }
0417     }
0418     return {trackIdx, quality};
0419   };
0420
0421   // Creating the map from TrackingParticle to SimTrackstersFromCP
0422   auto& simTrackstersFromCP = *result_fromCP;
0423   for (unsigned int i = 0; i < simTrackstersFromCP.size(); ++i) {
0424     if (simTrackstersFromCP[i].vertices().empty())
0425       continue;
0426     const auto& simTrack = caloparticles[simTrackstersFromCP[i].seedIndex()].g4Tracks()[0];
0427     UniqueSimTrackId simTkIds(simTrack.trackId(), simTrack.eventId());
0428     auto bestAssociatedRecoTrack = simTrackToRecoTrack(simTkIds);
0429     if (bestAssociatedRecoTrack.first != -1 and bestAssociatedRecoTrack.second > qualityCutTrack_) {
0430       auto trackIndex = bestAssociatedRecoTrack.first;
0431       simTrackstersFromCP[i].setTrackIdx(trackIndex);
0432     }
0433   }
0434
0435   auto& simTracksters = *result;
0436   // Creating the map from TrackingParticle to SimTrackster
0437   std::unordered_map<unsigned int, std::vector<unsigned int>> TPtoSimTracksterMap;
0438   for (unsigned int i = 0; i < simTracksters.size(); ++i) {
0439     const auto& simTrack = (simTracksters[i].seedID() == caloParticles_h.id())
0440                                ? caloparticles[simTracksters[i].seedIndex()].g4Tracks()[0]
0441                                : simclusters[simTracksters[i].seedIndex()].g4Tracks()[0];
0442     UniqueSimTrackId simTkIds(simTrack.trackId(), simTrack.eventId());
0443     auto bestAssociatedRecoTrack = simTrackToRecoTrack(simTkIds);
0444     if (bestAssociatedRecoTrack.first != -1 and bestAssociatedRecoTrack.second > qualityCutTrack_) {
0445       auto trackIndex = bestAssociatedRecoTrack.first;
0446       simTracksters[i].setTrackIdx(trackIndex);
0447     }
0448   }
0449
0450   edm::OrphanHandle<std::vector<Trackster>> simTracksters_h = evt.put(std::move(result));
0451
0452   // map between simTrack and Mtd SimTracksters to loop on them only one
0453   std::unordered_map<unsigned int, const MtdSimTrackster*> SimTrackToMtdST;
0454   for (unsigned int i = 0; i < MTDSimTracksters_h->size(); ++i) {
0455     const auto& simTrack = (*MTDSimTracksters_h)[i].g4Tracks()[0];
0456     SimTrackToMtdST[simTrack.trackId()] = &((*MTDSimTracksters_h)[i]);
0457   }
0458
0459   result_ticlCandidates->resize(result_fromCP->size());
0460   std::vector<int> toKeep;
0461   for (size_t i = 0; i < simTracksters_h->size(); ++i) {
0462     const auto& simTrackster = (*simTracksters_h)[i];
0463     int cp_index = (simTrackster.seedID() == caloParticles_h.id())
0464                        ? simTrackster.seedIndex()
0465                        : SimClusterToCaloParticleMap[simTrackster.seedIndex()];
0466     auto const& tCP = (*result_fromCP)[cp_index];
0467     if (!tCP.vertices().empty()) {
0468       auto trackIndex = tCP.trackIdx();
0469
0470       auto& cand = (*result_ticlCandidates)[cp_index];
0471       cand.addTrackster(edm::Ptr<Trackster>(simTracksters_h, i));
0472       if (trackIndex != -1 && caloparticles[cp_index].charge() != 0)
0473         cand.setTrackPtr(edm::Ptr<reco::Track>(recoTracks_h, trackIndex));
0474       toKeep.push_back(cp_index);
0475     }
0476   }
0477
0478   auto isHad = [](int pdgId) {
0479     pdgId = std::abs(pdgId);
0480     if (pdgId == 111)
0481       return false;
0482     return (pdgId > 100 and pdgId < 900) or (pdgId > 1000 and pdgId < 9000);
0483   };
0484
0485   for (size_t i = 0; i < result_ticlCandidates->size(); ++i) {
0486     auto cp_index = (*result_fromCP)[i].seedIndex();
0487     if (cp_index < 0)
0488       continue;
0489     auto& cand = (*result_ticlCandidates)[i];
0490     const auto& cp = caloparticles[cp_index];
0491     float rawEnergy = 0.f;
0492     float regressedEnergy = 0.f;
0493
0494     const auto& simTrack = cp.g4Tracks()[0];
0495     auto pos = SimTrackToMtdST.find(simTrack.trackId());
0496     if (pos != SimTrackToMtdST.end()) {
0497       auto MTDst = pos->second;
0498       // TODO: once the associators have been implemented check if the MTDst is associated with a reco before adding the MTD time
0499       cand.setMTDTime(MTDst->time(), 0);
0500     }
0501
0502     cand.setTime(cp.simTime(), 0);
0503
0504     for (const auto& trackster : cand.tracksters()) {
0505       rawEnergy += trackster->raw_energy();
0506       regressedEnergy += trackster->regressed_energy();
0507     }
0508     cand.setRawEnergy(rawEnergy);
0509
0510     auto pdgId = cp.pdgId();
0511     auto charge = cp.charge();
0512     if (cand.trackPtr().isNonnull()) {
0513       auto const& track = cand.trackPtr().get();
0514       if (std::abs(pdgId) == 13) {
0515         cand.setPdgId(pdgId);
0516       } else {
0517         cand.setPdgId((isHad(pdgId) ? 211 : 11) * charge);
0518       }
0519       cand.setCharge(charge);
0520       math::XYZTLorentzVector p4(regressedEnergy * track->momentum().unit().x(),
0521                                  regressedEnergy * track->momentum().unit().y(),
0522                                  regressedEnergy * track->momentum().unit().z(),
0523                                  regressedEnergy);
0524       cand.setP4(p4);
0525     } else {  // neutral candidates
0526       // a neutral candidate with a charged CaloParticle is charged without a reco track associated with it
0527       // set the charge = 0, but keep the real pdgId to keep track of that
0528       if (charge != 0)
0529         cand.setPdgId(isHad(pdgId) ? 211 : 11);
0530       else if (pdgId == 111)
0531         cand.setPdgId(pdgId);
0532       else
0533         cand.setPdgId(isHad(pdgId) ? 130 : 22);
0534       cand.setCharge(0);
0535
0536       auto particleType = tracksterParticleTypeFromPdgId(cand.pdgId(), 1);
0537       cand.setIdProbability(particleType, 1.f);
0538
0539       const auto& simTracksterFromCP = (*result_fromCP)[i];
0540       float regressedEnergy = simTracksterFromCP.regressed_energy();
0541       math::XYZTLorentzVector p4(regressedEnergy * simTracksterFromCP.barycenter().unit().x(),
0542                                  regressedEnergy * simTracksterFromCP.barycenter().unit().y(),
0543                                  regressedEnergy * simTracksterFromCP.barycenter().unit().z(),
0544                                  regressedEnergy);
0545       cand.setP4(p4);
0546     }
0547   }
0548
0549   std::vector<int> all_nums(result_fromCP->size());  // vector containing all caloparticles indexes
0550   std::iota(all_nums.begin(), all_nums.end(), 0);    // fill the vector with consecutive numbers starting from 0
0551
0552   std::vector<int> toRemove;
0553   std::set_difference(all_nums.begin(), all_nums.end(), toKeep.begin(), toKeep.end(), std::back_inserter(toRemove));
0554   std::sort(toRemove.begin(), toRemove.end(), [](int x, int y) { return x > y; });
0555   for (auto const& r : toRemove) {
0556     result_fromCP->erase(result_fromCP->begin() + r);
0557     result_ticlCandidates->erase(result_ticlCandidates->begin() + r);
0558   }
0559   evt.put(std::move(result_ticlCandidates));
0560   evt.put(std::move(output_mask));
0561   evt.put(std::move(result_fromCP), "fromCPs");
0562   evt.put(std::move(resultPU), "PU");
0563   evt.put(std::move(output_mask_fromCP), "fromCPs");
0564   evt.put(std::move(cpToSc_SimTrackstersMap));
0565 }