Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-29 11:58:06

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 "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
0032 
0033 #include "SimDataFormats/Vertex/interface/SimVertex.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 "RecoHGCal/TICL/interface/commons.h"
0041 
0042 #include "TrackstersPCA.h"
0043 #include <vector>
0044 #include <map>
0045 #include <iterator>
0046 #include <algorithm>
0047 #include <numeric>
0048 
0049 using namespace ticl;
0050 
0051 class SimTrackstersProducer : public edm::stream::EDProducer<> {
0052 public:
0053   explicit SimTrackstersProducer(const edm::ParameterSet&);
0054   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0055 
0056   void produce(edm::Event&, const edm::EventSetup&) override;
0057   void makePUTrackster(const std::vector<float>& inputClusterMask,
0058                        std::vector<float>& output_mask,
0059                        std::vector<Trackster>& result,
0060                        const edm::ProductID seed,
0061                        int loop_index);
0062 
0063   void addTrackster(const int index,
0064                     const std::vector<std::pair<edm::Ref<reco::CaloClusterCollection>, std::pair<float, float>>>& lcVec,
0065                     const std::vector<float>& inputClusterMask,
0066                     const float fractionCut_,
0067                     const float energy,
0068                     const int pdgId,
0069                     const int charge,
0070                     float time,
0071                     const edm::ProductID seed,
0072                     const Trackster::IterationIndex iter,
0073                     std::vector<float>& output_mask,
0074                     std::vector<Trackster>& result,
0075                     int& loop_index,
0076                     const bool add = false);
0077 
0078 private:
0079   std::string detector_;
0080   const bool doNose_ = false;
0081   const edm::EDGetTokenT<std::vector<reco::CaloCluster>> clusters_token_;
0082   const edm::EDGetTokenT<edm::ValueMap<std::pair<float, float>>> clustersTime_token_;
0083   const edm::EDGetTokenT<std::vector<float>> filtered_layerclusters_mask_token_;
0084 
0085   const edm::EDGetTokenT<std::vector<SimCluster>> simclusters_token_;
0086   const edm::EDGetTokenT<std::vector<CaloParticle>> caloparticles_token_;
0087 
0088   const edm::EDGetTokenT<hgcal::SimToRecoCollectionWithSimClusters> associatorMapSimClusterToReco_token_;
0089   const edm::EDGetTokenT<hgcal::SimToRecoCollection> associatorMapCaloParticleToReco_token_;
0090   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geom_token_;
0091   hgcal::RecHitTools rhtools_;
0092   const float fractionCut_;
0093   const float qualityCutTrack_;
0094   const edm::EDGetTokenT<std::vector<TrackingParticle>> trackingParticleToken_;
0095   const edm::EDGetTokenT<std::vector<SimVertex>> simVerticesToken_;
0096 
0097   const edm::EDGetTokenT<std::vector<reco::Track>> recoTracksToken_;
0098   const StringCutObjectSelector<reco::Track> cutTk_;
0099 
0100   const edm::EDGetTokenT<reco::SimToRecoCollection> associatormapStRsToken_;
0101   const edm::EDGetTokenT<reco::RecoToSimCollection> associatormapRtSsToken_;
0102   const edm::EDGetTokenT<SimTrackToTPMap> associationSimTrackToTPToken_;
0103 };
0104 
0105 DEFINE_FWK_MODULE(SimTrackstersProducer);
0106 
0107 SimTrackstersProducer::SimTrackstersProducer(const edm::ParameterSet& ps)
0108     : detector_(ps.getParameter<std::string>("detector")),
0109       doNose_(detector_ == "HFNose"),
0110       clusters_token_(consumes(ps.getParameter<edm::InputTag>("layer_clusters"))),
0111       clustersTime_token_(consumes(ps.getParameter<edm::InputTag>("time_layerclusters"))),
0112       filtered_layerclusters_mask_token_(consumes(ps.getParameter<edm::InputTag>("filtered_mask"))),
0113       simclusters_token_(consumes(ps.getParameter<edm::InputTag>("simclusters"))),
0114       caloparticles_token_(consumes(ps.getParameter<edm::InputTag>("caloparticles"))),
0115       associatorMapSimClusterToReco_token_(
0116           consumes(ps.getParameter<edm::InputTag>("layerClusterSimClusterAssociator"))),
0117       associatorMapCaloParticleToReco_token_(
0118           consumes(ps.getParameter<edm::InputTag>("layerClusterCaloParticleAssociator"))),
0119       geom_token_(esConsumes()),
0120       fractionCut_(ps.getParameter<double>("fractionCut")),
0121       qualityCutTrack_(ps.getParameter<double>("qualityCutTrack")),
0122       trackingParticleToken_(
0123           consumes<std::vector<TrackingParticle>>(ps.getParameter<edm::InputTag>("trackingParticles"))),
0124       simVerticesToken_(consumes<std::vector<SimVertex>>(ps.getParameter<edm::InputTag>("simVertices"))),
0125       recoTracksToken_(consumes<std::vector<reco::Track>>(ps.getParameter<edm::InputTag>("recoTracks"))),
0126       cutTk_(ps.getParameter<std::string>("cutTk")),
0127       associatormapStRsToken_(consumes(ps.getParameter<edm::InputTag>("tpToTrack"))),
0128       associationSimTrackToTPToken_(consumes(ps.getParameter<edm::InputTag>("simTrackToTPMap"))) {
0129   produces<TracksterCollection>();
0130   produces<std::vector<float>>();
0131   produces<TracksterCollection>("fromCPs");
0132   produces<TracksterCollection>("PU");
0133   produces<std::vector<float>>("fromCPs");
0134   produces<std::map<uint, std::vector<uint>>>();
0135   produces<std::vector<TICLCandidate>>();
0136 }
0137 
0138 void SimTrackstersProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0139   edm::ParameterSetDescription desc;
0140   desc.add<std::string>("detector", "HGCAL");
0141   desc.add<edm::InputTag>("layer_clusters", edm::InputTag("hgcalMergeLayerClusters"));
0142   desc.add<edm::InputTag>("time_layerclusters", edm::InputTag("hgcalMergeLayerClusters", "timeLayerCluster"));
0143   desc.add<edm::InputTag>("filtered_mask", edm::InputTag("filteredLayerClustersSimTracksters", "ticlSimTracksters"));
0144   desc.add<edm::InputTag>("simclusters", edm::InputTag("mix", "MergedCaloTruth"));
0145   desc.add<edm::InputTag>("caloparticles", edm::InputTag("mix", "MergedCaloTruth"));
0146   desc.add<edm::InputTag>("layerClusterSimClusterAssociator",
0147                           edm::InputTag("layerClusterSimClusterAssociationProducer"));
0148   desc.add<edm::InputTag>("layerClusterCaloParticleAssociator",
0149                           edm::InputTag("layerClusterCaloParticleAssociationProducer"));
0150   desc.add<edm::InputTag>("recoTracks", edm::InputTag("generalTracks"));
0151   desc.add<std::string>("cutTk",
0152                         "1.48 < abs(eta) < 3.0 && pt > 1. && quality(\"highPurity\") && "
0153                         "hitPattern().numberOfLostHits(\"MISSING_OUTER_HITS\") < 5");
0154   desc.add<edm::InputTag>("tpToTrack", edm::InputTag("trackingParticleRecoTrackAsssociation"));
0155 
0156   desc.add<edm::InputTag>("trackingParticles", edm::InputTag("mix", "MergedTrackTruth"));
0157   desc.add<edm::InputTag>("simVertices", edm::InputTag("g4SimHits"));
0158 
0159   desc.add<edm::InputTag>("simTrackToTPMap", edm::InputTag("simHitTPAssocProducer", "simTrackToTP"));
0160   desc.add<double>("fractionCut", 0.);
0161   desc.add<double>("qualityCutTrack", 0.75);
0162 
0163   descriptions.addWithDefaultLabel(desc);
0164 }
0165 void SimTrackstersProducer::makePUTrackster(const std::vector<float>& inputClusterMask,
0166                                             std::vector<float>& output_mask,
0167                                             std::vector<Trackster>& result,
0168                                             const edm::ProductID seed,
0169                                             int loop_index) {
0170   Trackster tmpTrackster;
0171   for (size_t i = 0; i < output_mask.size(); i++) {
0172     const float remaining_fraction = output_mask[i];
0173     if (remaining_fraction > 0.) {
0174       tmpTrackster.vertices().push_back(i);
0175       tmpTrackster.vertex_multiplicity().push_back(1. / remaining_fraction);
0176     }
0177   }
0178   tmpTrackster.setSeed(seed, 0);
0179   result.push_back(tmpTrackster);
0180 }
0181 
0182 void SimTrackstersProducer::addTrackster(
0183     const int index,
0184     const std::vector<std::pair<edm::Ref<reco::CaloClusterCollection>, std::pair<float, float>>>& lcVec,
0185     const std::vector<float>& inputClusterMask,
0186     const float fractionCut_,
0187     const float energy,
0188     const int pdgId,
0189     const int charge,
0190     const float time,
0191     const edm::ProductID seed,
0192     const Trackster::IterationIndex iter,
0193     std::vector<float>& output_mask,
0194     std::vector<Trackster>& result,
0195     int& loop_index,
0196     const bool add) {
0197   Trackster tmpTrackster;
0198   if (lcVec.empty()) {
0199     result[index] = tmpTrackster;
0200     return;
0201   }
0202 
0203   tmpTrackster.vertices().reserve(lcVec.size());
0204   tmpTrackster.vertex_multiplicity().reserve(lcVec.size());
0205   for (auto const& [lc, energyScorePair] : lcVec) {
0206     if (inputClusterMask[lc.index()] > 0) {
0207       double fraction = energyScorePair.first / lc->energy();
0208       if (fraction < fractionCut_)
0209         continue;
0210       tmpTrackster.vertices().push_back(lc.index());
0211       output_mask[lc.index()] -= fraction;
0212       tmpTrackster.vertex_multiplicity().push_back(1. / fraction);
0213     }
0214   }
0215 
0216   tmpTrackster.setIdProbability(tracksterParticleTypeFromPdgId(pdgId, charge), 1.f);
0217   tmpTrackster.setRegressedEnergy(energy);
0218   tmpTrackster.setIteration(iter);
0219   tmpTrackster.setSeed(seed, index);
0220   if (add) {
0221     result[index] = tmpTrackster;
0222     loop_index += 1;
0223   } else {
0224     result.push_back(tmpTrackster);
0225   }
0226 }
0227 
0228 void SimTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0229   auto result = std::make_unique<TracksterCollection>();
0230   auto output_mask = std::make_unique<std::vector<float>>();
0231   auto result_fromCP = std::make_unique<TracksterCollection>();
0232   auto resultPU = std::make_unique<TracksterCollection>();
0233   auto output_mask_fromCP = std::make_unique<std::vector<float>>();
0234   auto cpToSc_SimTrackstersMap = std::make_unique<std::map<uint, std::vector<uint>>>();
0235   auto result_ticlCandidates = std::make_unique<std::vector<TICLCandidate>>();
0236 
0237   const auto& layerClusters = evt.get(clusters_token_);
0238   const auto& layerClustersTimes = evt.get(clustersTime_token_);
0239   const auto& inputClusterMask = evt.get(filtered_layerclusters_mask_token_);
0240   output_mask->resize(layerClusters.size(), 1.f);
0241   output_mask_fromCP->resize(layerClusters.size(), 1.f);
0242 
0243   const auto& simclusters = evt.get(simclusters_token_);
0244   edm::Handle<std::vector<CaloParticle>> caloParticles_h;
0245   evt.getByToken(caloparticles_token_, caloParticles_h);
0246   const auto& caloparticles = *caloParticles_h;
0247 
0248   const auto& simClustersToRecoColl = evt.get(associatorMapSimClusterToReco_token_);
0249   const auto& caloParticlesToRecoColl = evt.get(associatorMapCaloParticleToReco_token_);
0250   const auto& simVertices = evt.get(simVerticesToken_);
0251 
0252   edm::Handle<std::vector<TrackingParticle>> trackingParticles_h;
0253   evt.getByToken(trackingParticleToken_, trackingParticles_h);
0254   edm::Handle<std::vector<reco::Track>> recoTracks_h;
0255   evt.getByToken(recoTracksToken_, recoTracks_h);
0256   const auto& TPtoRecoTrackMap = evt.get(associatormapStRsToken_);
0257   const auto& simTrackToTPMap = evt.get(associationSimTrackToTPToken_);
0258   const auto& recoTracks = *recoTracks_h;
0259 
0260   const auto& geom = es.getData(geom_token_);
0261   rhtools_.setGeometry(geom);
0262   const auto num_simclusters = simclusters.size();
0263   result->reserve(num_simclusters);  // Conservative size, will call shrink_to_fit later
0264   const auto num_caloparticles = caloparticles.size();
0265   result_fromCP->resize(num_caloparticles);
0266   std::map<uint, uint> SimClusterToCaloParticleMap;
0267   int loop_index = 0;
0268   for (const auto& [key, lcVec] : caloParticlesToRecoColl) {
0269     auto const& cp = *(key);
0270     auto cpIndex = &cp - &caloparticles[0];
0271     for (const auto& scRef : cp.simClusters()) {
0272       auto const& sc = *(scRef);
0273       auto const scIndex = &sc - &simclusters[0];
0274       SimClusterToCaloParticleMap[scIndex] = cpIndex;
0275     }
0276 
0277     auto regr_energy = cp.energy();
0278     std::vector<uint> scSimTracksterIdx;
0279     scSimTracksterIdx.reserve(cp.simClusters().size());
0280 
0281     // Create a Trackster from the object entering HGCal
0282     if (cp.g4Tracks()[0].crossedBoundary()) {
0283       regr_energy = cp.g4Tracks()[0].getMomentumAtBoundary().energy();
0284       float time = cp.g4Tracks()[0].getPositionAtBoundary().t();
0285       addTrackster(cpIndex,
0286                    lcVec,
0287                    inputClusterMask,
0288                    fractionCut_,
0289                    regr_energy,
0290                    cp.pdgId(),
0291                    cp.charge(),
0292                    time,
0293                    key.id(),
0294                    ticl::Trackster::SIM,
0295                    *output_mask,
0296                    *result,
0297                    loop_index);
0298     } else {
0299       for (const auto& scRef : cp.simClusters()) {
0300         const auto& it = simClustersToRecoColl.find(scRef);
0301         if (it == simClustersToRecoColl.end())
0302           continue;
0303         const auto& lcVec = it->val;
0304         auto const& sc = *(scRef);
0305         auto const scIndex = &sc - &simclusters[0];
0306 
0307         addTrackster(scIndex,
0308                      lcVec,
0309                      inputClusterMask,
0310                      fractionCut_,
0311                      sc.g4Tracks()[0].getMomentumAtBoundary().energy(),
0312                      sc.pdgId(),
0313                      sc.charge(),
0314                      sc.g4Tracks()[0].getPositionAtBoundary().t(),
0315                      scRef.id(),
0316                      ticl::Trackster::SIM,
0317                      *output_mask,
0318                      *result,
0319                      loop_index);
0320 
0321         if (result->empty())
0322           continue;
0323         const auto index = result->size() - 1;
0324         if (std::find(scSimTracksterIdx.begin(), scSimTracksterIdx.end(), index) == scSimTracksterIdx.end()) {
0325           scSimTracksterIdx.emplace_back(index);
0326         }
0327       }
0328       scSimTracksterIdx.shrink_to_fit();
0329     }
0330     float time = simVertices[cp.g4Tracks()[0].vertIndex()].position().t();
0331     // Create a Trackster from any CP
0332     addTrackster(cpIndex,
0333                  lcVec,
0334                  inputClusterMask,
0335                  fractionCut_,
0336                  regr_energy,
0337                  cp.pdgId(),
0338                  cp.charge(),
0339                  time,
0340                  key.id(),
0341                  ticl::Trackster::SIM_CP,
0342                  *output_mask_fromCP,
0343                  *result_fromCP,
0344                  loop_index,
0345                  true);
0346 
0347     if (result_fromCP->empty())
0348       continue;
0349     const auto index = loop_index - 1;
0350     if (cpToSc_SimTrackstersMap->find(index) == cpToSc_SimTrackstersMap->end()) {
0351       (*cpToSc_SimTrackstersMap)[index] = scSimTracksterIdx;
0352     }
0353   }
0354   // TODO: remove time computation from PCA calculation and
0355   //       store time from boundary position in simTracksters
0356   ticl::assignPCAtoTracksters(
0357       *result, layerClusters, layerClustersTimes, rhtools_.getPositionLayer(rhtools_.lastLayerEE(doNose_)).z());
0358   result->shrink_to_fit();
0359   ticl::assignPCAtoTracksters(
0360       *result_fromCP, layerClusters, layerClustersTimes, rhtools_.getPositionLayer(rhtools_.lastLayerEE(doNose_)).z());
0361 
0362   makePUTrackster(inputClusterMask, *output_mask, *resultPU, caloParticles_h.id(), 0);
0363 
0364   auto simTrackToRecoTrack = [&](UniqueSimTrackId simTkId) -> std::pair<int, float> {
0365     int trackIdx = -1;
0366     float quality = 0.f;
0367     auto ipos = simTrackToTPMap.mapping.find(simTkId);
0368     if (ipos != simTrackToTPMap.mapping.end()) {
0369       auto jpos = TPtoRecoTrackMap.find((ipos->second));
0370       if (jpos != TPtoRecoTrackMap.end()) {
0371         auto& associatedRecoTracks = jpos->val;
0372         if (!associatedRecoTracks.empty()) {
0373           // associated reco tracks are sorted by decreasing quality
0374           if (associatedRecoTracks[0].second > qualityCutTrack_) {
0375             trackIdx = &(*associatedRecoTracks[0].first) - &recoTracks[0];
0376             quality = associatedRecoTracks[0].second;
0377           }
0378         }
0379       }
0380     }
0381     return {trackIdx, quality};
0382   };
0383 
0384   // Creating the map from TrackingParticle to SimTrackstersFromCP
0385   auto& simTrackstersFromCP = *result_fromCP;
0386   for (unsigned int i = 0; i < simTrackstersFromCP.size(); ++i) {
0387     if (simTrackstersFromCP[i].vertices().empty())
0388       continue;
0389     const auto& simTrack = caloparticles[simTrackstersFromCP[i].seedIndex()].g4Tracks()[0];
0390     UniqueSimTrackId simTkIds(simTrack.trackId(), simTrack.eventId());
0391     auto bestAssociatedRecoTrack = simTrackToRecoTrack(simTkIds);
0392     if (bestAssociatedRecoTrack.first != -1 and bestAssociatedRecoTrack.second > qualityCutTrack_) {
0393       auto trackIndex = bestAssociatedRecoTrack.first;
0394       simTrackstersFromCP[i].setTrackIdx(trackIndex);
0395     }
0396   }
0397 
0398   auto& simTracksters = *result;
0399   // Creating the map from TrackingParticle to SimTrackster
0400   std::unordered_map<unsigned int, std::vector<unsigned int>> TPtoSimTracksterMap;
0401   for (unsigned int i = 0; i < simTracksters.size(); ++i) {
0402     const auto& simTrack = (simTracksters[i].seedID() == caloParticles_h.id())
0403                                ? caloparticles[simTracksters[i].seedIndex()].g4Tracks()[0]
0404                                : simclusters[simTracksters[i].seedIndex()].g4Tracks()[0];
0405     UniqueSimTrackId simTkIds(simTrack.trackId(), simTrack.eventId());
0406     auto bestAssociatedRecoTrack = simTrackToRecoTrack(simTkIds);
0407     if (bestAssociatedRecoTrack.first != -1 and bestAssociatedRecoTrack.second > qualityCutTrack_) {
0408       auto trackIndex = bestAssociatedRecoTrack.first;
0409       simTracksters[i].setTrackIdx(trackIndex);
0410     }
0411   }
0412 
0413   edm::OrphanHandle<std::vector<Trackster>> simTracksters_h = evt.put(std::move(result));
0414 
0415   result_ticlCandidates->resize(result_fromCP->size());
0416   std::vector<int> toKeep;
0417   for (size_t i = 0; i < simTracksters_h->size(); ++i) {
0418     const auto& simTrackster = (*simTracksters_h)[i];
0419     int cp_index = (simTrackster.seedID() == caloParticles_h.id())
0420                        ? simTrackster.seedIndex()
0421                        : SimClusterToCaloParticleMap[simTrackster.seedIndex()];
0422     auto const& tCP = (*result_fromCP)[cp_index];
0423     if (!tCP.vertices().empty()) {
0424       auto trackIndex = tCP.trackIdx();
0425 
0426       auto& cand = (*result_ticlCandidates)[cp_index];
0427       cand.addTrackster(edm::Ptr<Trackster>(simTracksters_h, i));
0428       if (trackIndex != -1 and (trackIndex < 0 or trackIndex >= (long int)recoTracks.size())) {
0429       }
0430       cand.setTime((*result_fromCP)[cp_index].time());
0431       cand.setTimeError(0);
0432       if (trackIndex != -1 && caloparticles[cp_index].charge() != 0)
0433         cand.setTrackPtr(edm::Ptr<reco::Track>(recoTracks_h, trackIndex));
0434       toKeep.push_back(cp_index);
0435     }
0436   }
0437 
0438   auto isHad = [](int pdgId) {
0439     pdgId = std::abs(pdgId);
0440     if (pdgId == 111)
0441       return false;
0442     return (pdgId > 100 and pdgId < 900) or (pdgId > 1000 and pdgId < 9000);
0443   };
0444 
0445   for (size_t i = 0; i < result_ticlCandidates->size(); ++i) {
0446     auto cp_index = (*result_fromCP)[i].seedIndex();
0447     if (cp_index < 0)
0448       continue;
0449     auto& cand = (*result_ticlCandidates)[i];
0450     const auto& cp = caloparticles[cp_index];
0451     float rawEnergy = 0.f;
0452     float regressedEnergy = 0.f;
0453 
0454     for (const auto& trackster : cand.tracksters()) {
0455       rawEnergy += trackster->raw_energy();
0456       regressedEnergy += trackster->regressed_energy();
0457     }
0458     cand.setRawEnergy(rawEnergy);
0459 
0460     auto pdgId = cp.pdgId();
0461     auto charge = cp.charge();
0462     if (cand.trackPtr().isNonnull() and charge == 0) {
0463     }
0464     if (cand.trackPtr().isNonnull() and charge != 0) {
0465       auto const& track = cand.trackPtr().get();
0466       if (std::abs(pdgId) == 13) {
0467         cand.setPdgId(pdgId);
0468       } else {
0469         cand.setPdgId((isHad(pdgId) ? 211 : 11) * charge);
0470       }
0471       cand.setCharge(charge);
0472       math::XYZTLorentzVector p4(regressedEnergy * track->momentum().unit().x(),
0473                                  regressedEnergy * track->momentum().unit().y(),
0474                                  regressedEnergy * track->momentum().unit().z(),
0475                                  regressedEnergy);
0476       cand.setP4(p4);
0477     } else {  // neutral candidates
0478       cand.setPdgId(isHad(pdgId) ? 130 : 22);
0479       cand.setCharge(0);
0480 
0481       auto particleType = tracksterParticleTypeFromPdgId(cand.pdgId(), 1);
0482       cand.setIdProbability(particleType, 1.f);
0483 
0484       const auto& simTracksterFromCP = (*result_fromCP)[i];
0485       float regressedEnergy = simTracksterFromCP.regressed_energy();
0486       math::XYZTLorentzVector p4(regressedEnergy * simTracksterFromCP.barycenter().unit().x(),
0487                                  regressedEnergy * simTracksterFromCP.barycenter().unit().y(),
0488                                  regressedEnergy * simTracksterFromCP.barycenter().unit().z(),
0489                                  regressedEnergy);
0490       cand.setP4(p4);
0491     }
0492   }
0493 
0494   std::vector<int> all_nums(result_fromCP->size());  // vector containing all caloparticles indexes
0495   std::iota(all_nums.begin(), all_nums.end(), 0);    // fill the vector with consecutive numbers starting from 0
0496 
0497   std::vector<int> toRemove;
0498   std::set_difference(all_nums.begin(), all_nums.end(), toKeep.begin(), toKeep.end(), std::back_inserter(toRemove));
0499   std::sort(toRemove.begin(), toRemove.end(), [](int x, int y) { return x > y; });
0500   for (auto const& r : toRemove) {
0501     result_fromCP->erase(result_fromCP->begin() + r);
0502     result_ticlCandidates->erase(result_ticlCandidates->begin() + r);
0503   }
0504   evt.put(std::move(result_ticlCandidates));
0505   evt.put(std::move(output_mask));
0506   evt.put(std::move(result_fromCP), "fromCPs");
0507   evt.put(std::move(resultPU), "PU");
0508   evt.put(std::move(output_mask_fromCP), "fromCPs");
0509   evt.put(std::move(cpToSc_SimTrackstersMap));
0510 }