File indexing completed on 2024-10-16 05:06:31
0001
0002
0003
0004
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/Vertex/interface/SimVertex.h"
0036
0037 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0038 #include "SimDataFormats/TrackingAnalysis/interface/UniqueSimTrackId.h"
0039
0040 #include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h"
0041
0042 #include "DataFormats/HGCalReco/interface/Common.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 const edm::EDGetTokenT<std::vector<SimVertex>> simVerticesToken_;
0100
0101 const edm::EDGetTokenT<std::vector<reco::Track>> recoTracksToken_;
0102 const StringCutObjectSelector<reco::Track> cutTk_;
0103
0104 const edm::EDGetTokenT<reco::SimToRecoCollection> associatormapStRsToken_;
0105 const edm::EDGetTokenT<reco::RecoToSimCollection> associatormapRtSsToken_;
0106 const edm::EDGetTokenT<SimTrackToTPMap> associationSimTrackToTPToken_;
0107 };
0108
0109 DEFINE_FWK_MODULE(SimTrackstersProducer);
0110
0111 SimTrackstersProducer::SimTrackstersProducer(const edm::ParameterSet& ps)
0112 : detector_(ps.getParameter<std::string>("detector")),
0113 doNose_(detector_ == "HFNose"),
0114 computeLocalTime_(ps.getParameter<bool>("computeLocalTime")),
0115 clusters_token_(consumes(ps.getParameter<edm::InputTag>("layer_clusters"))),
0116 clustersTime_token_(consumes(ps.getParameter<edm::InputTag>("time_layerclusters"))),
0117 filtered_layerclusters_mask_token_(consumes(ps.getParameter<edm::InputTag>("filtered_mask"))),
0118 simclusters_token_(consumes(ps.getParameter<edm::InputTag>("simclusters"))),
0119 caloparticles_token_(consumes(ps.getParameter<edm::InputTag>("caloparticles"))),
0120 MTDSimTrackstersToken_(consumes<MtdSimTracksterCollection>(ps.getParameter<edm::InputTag>("MtdSimTracksters"))),
0121 associatorMapSimClusterToReco_token_(
0122 consumes(ps.getParameter<edm::InputTag>("layerClusterSimClusterAssociator"))),
0123 associatorMapCaloParticleToReco_token_(
0124 consumes(ps.getParameter<edm::InputTag>("layerClusterCaloParticleAssociator"))),
0125 geom_token_(esConsumes()),
0126 fractionCut_(ps.getParameter<double>("fractionCut")),
0127 qualityCutTrack_(ps.getParameter<double>("qualityCutTrack")),
0128 trackingParticleToken_(
0129 consumes<std::vector<TrackingParticle>>(ps.getParameter<edm::InputTag>("trackingParticles"))),
0130 simVerticesToken_(consumes<std::vector<SimVertex>>(ps.getParameter<edm::InputTag>("simVertices"))),
0131 recoTracksToken_(consumes<std::vector<reco::Track>>(ps.getParameter<edm::InputTag>("recoTracks"))),
0132 cutTk_(ps.getParameter<std::string>("cutTk")),
0133 associatormapStRsToken_(consumes(ps.getParameter<edm::InputTag>("tpToTrack"))),
0134 associationSimTrackToTPToken_(consumes(ps.getParameter<edm::InputTag>("simTrackToTPMap"))) {
0135 produces<TracksterCollection>();
0136 produces<std::vector<float>>();
0137 produces<TracksterCollection>("fromCPs");
0138 produces<TracksterCollection>("PU");
0139 produces<std::vector<float>>("fromCPs");
0140 produces<std::map<uint, std::vector<uint>>>();
0141 produces<std::vector<TICLCandidate>>();
0142 }
0143
0144 void SimTrackstersProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0145 edm::ParameterSetDescription desc;
0146 desc.add<std::string>("detector", "HGCAL");
0147 desc.add<bool>("computeLocalTime", "false");
0148 desc.add<edm::InputTag>("layer_clusters", edm::InputTag("hgcalMergeLayerClusters"));
0149 desc.add<edm::InputTag>("time_layerclusters", edm::InputTag("hgcalMergeLayerClusters", "timeLayerCluster"));
0150 desc.add<edm::InputTag>("filtered_mask", edm::InputTag("filteredLayerClustersSimTracksters", "ticlSimTracksters"));
0151 desc.add<edm::InputTag>("simclusters", edm::InputTag("mix", "MergedCaloTruth"));
0152 desc.add<edm::InputTag>("caloparticles", edm::InputTag("mix", "MergedCaloTruth"));
0153 desc.add<edm::InputTag>("MtdSimTracksters", edm::InputTag("mix", "MergedMtdTruthST"));
0154 desc.add<edm::InputTag>("layerClusterSimClusterAssociator",
0155 edm::InputTag("layerClusterSimClusterAssociationProducer"));
0156 desc.add<edm::InputTag>("layerClusterCaloParticleAssociator",
0157 edm::InputTag("layerClusterCaloParticleAssociationProducer"));
0158 desc.add<edm::InputTag>("recoTracks", edm::InputTag("generalTracks"));
0159 desc.add<std::string>("cutTk",
0160 "1.48 < abs(eta) < 3.0 && pt > 1. && quality(\"highPurity\") && "
0161 "hitPattern().numberOfLostHits(\"MISSING_OUTER_HITS\") < 5");
0162 desc.add<edm::InputTag>("tpToTrack", edm::InputTag("trackingParticleRecoTrackAsssociation"));
0163
0164 desc.add<edm::InputTag>("trackingParticles", edm::InputTag("mix", "MergedTrackTruth"));
0165 desc.add<edm::InputTag>("simVertices", edm::InputTag("g4SimHits"));
0166
0167 desc.add<edm::InputTag>("simTrackToTPMap", edm::InputTag("simHitTPAssocProducer", "simTrackToTP"));
0168 desc.add<double>("fractionCut", 0.);
0169 desc.add<double>("qualityCutTrack", 0.75);
0170
0171 descriptions.addWithDefaultLabel(desc);
0172 }
0173 void SimTrackstersProducer::makePUTrackster(const std::vector<float>& inputClusterMask,
0174 std::vector<float>& output_mask,
0175 std::vector<Trackster>& result,
0176 const edm::ProductID seed,
0177 int loop_index) {
0178 Trackster tmpTrackster;
0179 for (size_t i = 0; i < output_mask.size(); i++) {
0180 const float remaining_fraction = output_mask[i];
0181 if (remaining_fraction > std::numeric_limits<float>::epsilon()) {
0182 tmpTrackster.vertices().push_back(i);
0183 tmpTrackster.vertex_multiplicity().push_back(1. / remaining_fraction);
0184 }
0185 }
0186 tmpTrackster.setSeed(seed, 0);
0187 result.push_back(tmpTrackster);
0188 }
0189
0190 void SimTrackstersProducer::addTrackster(
0191 const int index,
0192 const std::vector<std::pair<edm::Ref<reco::CaloClusterCollection>, std::pair<float, float>>>& lcVec,
0193 const std::vector<float>& inputClusterMask,
0194 const float fractionCut_,
0195 const float energy,
0196 const int pdgId,
0197 const int charge,
0198 const float time,
0199 const edm::ProductID seed,
0200 const Trackster::IterationIndex iter,
0201 std::vector<float>& output_mask,
0202 std::vector<Trackster>& result,
0203 int& loop_index,
0204 const bool add) {
0205 Trackster tmpTrackster;
0206 if (lcVec.empty()) {
0207 result[index] = tmpTrackster;
0208 return;
0209 }
0210
0211 tmpTrackster.vertices().reserve(lcVec.size());
0212 tmpTrackster.vertex_multiplicity().reserve(lcVec.size());
0213 for (auto const& [lc, energyScorePair] : lcVec) {
0214 if (inputClusterMask[lc.index()] > 0) {
0215 float fraction = energyScorePair.first / lc->energy();
0216 if (fraction < fractionCut_)
0217 continue;
0218 tmpTrackster.vertices().push_back(lc.index());
0219 output_mask[lc.index()] -= fraction;
0220 tmpTrackster.vertex_multiplicity().push_back(1. / fraction);
0221 }
0222 }
0223
0224 tmpTrackster.setIdProbability(tracksterParticleTypeFromPdgId(pdgId, charge), 1.f);
0225 tmpTrackster.setRegressedEnergy(energy);
0226 tmpTrackster.setIteration(iter);
0227 tmpTrackster.setSeed(seed, index);
0228 tmpTrackster.setBoundaryTime(time * 1e9);
0229 if (add) {
0230 result[index] = tmpTrackster;
0231 loop_index += 1;
0232 } else {
0233 result.push_back(tmpTrackster);
0234 }
0235 }
0236
0237 void SimTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0238 auto result = std::make_unique<TracksterCollection>();
0239 auto output_mask = std::make_unique<std::vector<float>>();
0240 auto result_fromCP = std::make_unique<TracksterCollection>();
0241 auto resultPU = std::make_unique<TracksterCollection>();
0242 auto output_mask_fromCP = std::make_unique<std::vector<float>>();
0243 auto cpToSc_SimTrackstersMap = std::make_unique<std::map<uint, std::vector<uint>>>();
0244 auto result_ticlCandidates = std::make_unique<std::vector<TICLCandidate>>();
0245
0246 const auto& layerClusters = evt.get(clusters_token_);
0247 const auto& layerClustersTimes = evt.get(clustersTime_token_);
0248 const auto& inputClusterMask = evt.get(filtered_layerclusters_mask_token_);
0249 output_mask->resize(layerClusters.size(), 1.f);
0250 output_mask_fromCP->resize(layerClusters.size(), 1.f);
0251
0252 const auto& simclusters = evt.get(simclusters_token_);
0253 edm::Handle<std::vector<CaloParticle>> caloParticles_h;
0254 evt.getByToken(caloparticles_token_, caloParticles_h);
0255 const auto& caloparticles = *caloParticles_h;
0256
0257 edm::Handle<MtdSimTracksterCollection> MTDSimTracksters_h;
0258 evt.getByToken(MTDSimTrackstersToken_, MTDSimTracksters_h);
0259
0260 const auto& simClustersToRecoColl = evt.get(associatorMapSimClusterToReco_token_);
0261 const auto& caloParticlesToRecoColl = evt.get(associatorMapCaloParticleToReco_token_);
0262 const auto& simVertices = evt.get(simVerticesToken_);
0263
0264 edm::Handle<std::vector<TrackingParticle>> trackingParticles_h;
0265 evt.getByToken(trackingParticleToken_, trackingParticles_h);
0266 edm::Handle<std::vector<reco::Track>> recoTracks_h;
0267 evt.getByToken(recoTracksToken_, recoTracks_h);
0268 const auto& TPtoRecoTrackMap = evt.get(associatormapStRsToken_);
0269 const auto& simTrackToTPMap = evt.get(associationSimTrackToTPToken_);
0270 const auto& recoTracks = *recoTracks_h;
0271
0272 const auto& geom = es.getData(geom_token_);
0273 rhtools_.setGeometry(geom);
0274 const auto num_simclusters = simclusters.size();
0275 result->reserve(num_simclusters);
0276 const auto num_caloparticles = caloparticles.size();
0277 result_fromCP->resize(num_caloparticles);
0278 std::map<uint, uint> SimClusterToCaloParticleMap;
0279 int loop_index = 0;
0280 for (const auto& [key, lcVec] : caloParticlesToRecoColl) {
0281 auto const& cp = *(key);
0282 auto cpIndex = &cp - &caloparticles[0];
0283 for (const auto& scRef : cp.simClusters()) {
0284 auto const& sc = *(scRef);
0285 auto const scIndex = &sc - &simclusters[0];
0286 SimClusterToCaloParticleMap[scIndex] = cpIndex;
0287 }
0288
0289 auto regr_energy = cp.energy();
0290 std::vector<uint> scSimTracksterIdx;
0291 scSimTracksterIdx.reserve(cp.simClusters().size());
0292
0293
0294 if (cp.g4Tracks()[0].crossedBoundary()) {
0295 regr_energy = cp.g4Tracks()[0].getMomentumAtBoundary().energy();
0296 float time = cp.g4Tracks()[0].getPositionAtBoundary().t();
0297 addTrackster(cpIndex,
0298 lcVec,
0299 inputClusterMask,
0300 fractionCut_,
0301 regr_energy,
0302 cp.pdgId(),
0303 cp.charge(),
0304 time,
0305 key.id(),
0306 ticl::Trackster::SIM,
0307 *output_mask,
0308 *result,
0309 loop_index);
0310 } else {
0311 for (const auto& scRef : cp.simClusters()) {
0312 const auto& it = simClustersToRecoColl.find(scRef);
0313 if (it == simClustersToRecoColl.end())
0314 continue;
0315 const auto& lcVec = it->val;
0316 auto const& sc = *(scRef);
0317 auto const scIndex = &sc - &simclusters[0];
0318
0319 addTrackster(scIndex,
0320 lcVec,
0321 inputClusterMask,
0322 fractionCut_,
0323 sc.g4Tracks()[0].getMomentumAtBoundary().energy(),
0324 sc.pdgId(),
0325 sc.charge(),
0326 sc.g4Tracks()[0].getPositionAtBoundary().t(),
0327 scRef.id(),
0328 ticl::Trackster::SIM,
0329 *output_mask,
0330 *result,
0331 loop_index);
0332
0333 if (result->empty())
0334 continue;
0335 const auto index = result->size() - 1;
0336 if (std::find(scSimTracksterIdx.begin(), scSimTracksterIdx.end(), index) == scSimTracksterIdx.end()) {
0337 scSimTracksterIdx.emplace_back(index);
0338 }
0339 }
0340 scSimTracksterIdx.shrink_to_fit();
0341 }
0342 float time = simVertices[cp.g4Tracks()[0].vertIndex()].position().t();
0343
0344 addTrackster(cpIndex,
0345 lcVec,
0346 inputClusterMask,
0347 fractionCut_,
0348 regr_energy,
0349 cp.pdgId(),
0350 cp.charge(),
0351 time,
0352 key.id(),
0353 ticl::Trackster::SIM_CP,
0354 *output_mask_fromCP,
0355 *result_fromCP,
0356 loop_index,
0357 true);
0358
0359 if (result_fromCP->empty())
0360 continue;
0361 const auto index = loop_index - 1;
0362 if (cpToSc_SimTrackstersMap->find(index) == cpToSc_SimTrackstersMap->end()) {
0363 (*cpToSc_SimTrackstersMap)[index] = scSimTracksterIdx;
0364 }
0365 }
0366
0367
0368 ticl::assignPCAtoTracksters(*result,
0369 layerClusters,
0370 layerClustersTimes,
0371 rhtools_.getPositionLayer(rhtools_.lastLayerEE(doNose_)).z(),
0372 rhtools_,
0373 computeLocalTime_);
0374 result->shrink_to_fit();
0375 ticl::assignPCAtoTracksters(*result_fromCP,
0376 layerClusters,
0377 layerClustersTimes,
0378 rhtools_.getPositionLayer(rhtools_.lastLayerEE(doNose_)).z(),
0379 rhtools_,
0380 computeLocalTime_);
0381
0382 makePUTrackster(inputClusterMask, *output_mask, *resultPU, caloParticles_h.id(), 0);
0383
0384 auto simTrackToRecoTrack = [&](UniqueSimTrackId simTkId) -> std::pair<int, float> {
0385 int trackIdx = -1;
0386 float quality = 0.f;
0387 auto ipos = simTrackToTPMap.mapping.find(simTkId);
0388 if (ipos != simTrackToTPMap.mapping.end()) {
0389 auto jpos = TPtoRecoTrackMap.find((ipos->second));
0390 if (jpos != TPtoRecoTrackMap.end()) {
0391 auto& associatedRecoTracks = jpos->val;
0392 if (!associatedRecoTracks.empty()) {
0393
0394 if (associatedRecoTracks[0].second > qualityCutTrack_) {
0395 trackIdx = &(*associatedRecoTracks[0].first) - &recoTracks[0];
0396 quality = associatedRecoTracks[0].second;
0397 }
0398 }
0399 }
0400 }
0401 return {trackIdx, quality};
0402 };
0403
0404
0405 auto& simTrackstersFromCP = *result_fromCP;
0406 for (unsigned int i = 0; i < simTrackstersFromCP.size(); ++i) {
0407 if (simTrackstersFromCP[i].vertices().empty())
0408 continue;
0409 const auto& simTrack = caloparticles[simTrackstersFromCP[i].seedIndex()].g4Tracks()[0];
0410 UniqueSimTrackId simTkIds(simTrack.trackId(), simTrack.eventId());
0411 auto bestAssociatedRecoTrack = simTrackToRecoTrack(simTkIds);
0412 if (bestAssociatedRecoTrack.first != -1 and bestAssociatedRecoTrack.second > qualityCutTrack_) {
0413 auto trackIndex = bestAssociatedRecoTrack.first;
0414 simTrackstersFromCP[i].setTrackIdx(trackIndex);
0415 }
0416 }
0417
0418 auto& simTracksters = *result;
0419
0420 std::unordered_map<unsigned int, std::vector<unsigned int>> TPtoSimTracksterMap;
0421 for (unsigned int i = 0; i < simTracksters.size(); ++i) {
0422 const auto& simTrack = (simTracksters[i].seedID() == caloParticles_h.id())
0423 ? caloparticles[simTracksters[i].seedIndex()].g4Tracks()[0]
0424 : simclusters[simTracksters[i].seedIndex()].g4Tracks()[0];
0425 UniqueSimTrackId simTkIds(simTrack.trackId(), simTrack.eventId());
0426 auto bestAssociatedRecoTrack = simTrackToRecoTrack(simTkIds);
0427 if (bestAssociatedRecoTrack.first != -1 and bestAssociatedRecoTrack.second > qualityCutTrack_) {
0428 auto trackIndex = bestAssociatedRecoTrack.first;
0429 simTracksters[i].setTrackIdx(trackIndex);
0430 }
0431 }
0432
0433 edm::OrphanHandle<std::vector<Trackster>> simTracksters_h = evt.put(std::move(result));
0434
0435
0436 std::unordered_map<unsigned int, const MtdSimTrackster*> SimTrackToMtdST;
0437 for (unsigned int i = 0; i < MTDSimTracksters_h->size(); ++i) {
0438 const auto& simTrack = (*MTDSimTracksters_h)[i].g4Tracks()[0];
0439 SimTrackToMtdST[simTrack.trackId()] = &((*MTDSimTracksters_h)[i]);
0440 }
0441
0442 result_ticlCandidates->resize(result_fromCP->size());
0443 std::vector<int> toKeep;
0444 for (size_t i = 0; i < simTracksters_h->size(); ++i) {
0445 const auto& simTrackster = (*simTracksters_h)[i];
0446 int cp_index = (simTrackster.seedID() == caloParticles_h.id())
0447 ? simTrackster.seedIndex()
0448 : SimClusterToCaloParticleMap[simTrackster.seedIndex()];
0449 auto const& tCP = (*result_fromCP)[cp_index];
0450 if (!tCP.vertices().empty()) {
0451 auto trackIndex = tCP.trackIdx();
0452
0453 auto& cand = (*result_ticlCandidates)[cp_index];
0454 cand.addTrackster(edm::Ptr<Trackster>(simTracksters_h, i));
0455 if (trackIndex != -1 && caloparticles[cp_index].charge() != 0)
0456 cand.setTrackPtr(edm::Ptr<reco::Track>(recoTracks_h, trackIndex));
0457 toKeep.push_back(cp_index);
0458 }
0459 }
0460
0461 auto isHad = [](int pdgId) {
0462 pdgId = std::abs(pdgId);
0463 if (pdgId == 111)
0464 return false;
0465 return (pdgId > 100 and pdgId < 900) or (pdgId > 1000 and pdgId < 9000);
0466 };
0467
0468 for (size_t i = 0; i < result_ticlCandidates->size(); ++i) {
0469 auto cp_index = (*result_fromCP)[i].seedIndex();
0470 if (cp_index < 0)
0471 continue;
0472 auto& cand = (*result_ticlCandidates)[i];
0473 const auto& cp = caloparticles[cp_index];
0474 float rawEnergy = 0.f;
0475 float regressedEnergy = 0.f;
0476
0477 const auto& simTrack = cp.g4Tracks()[0];
0478 auto pos = SimTrackToMtdST.find(simTrack.trackId());
0479 if (pos != SimTrackToMtdST.end()) {
0480 auto MTDst = pos->second;
0481
0482 cand.setMTDTime(MTDst->time(), 0);
0483 }
0484
0485 cand.setTime(simVertices[cp.g4Tracks()[0].vertIndex()].position().t() * pow(10, 9), 0);
0486
0487 for (const auto& trackster : cand.tracksters()) {
0488 rawEnergy += trackster->raw_energy();
0489 regressedEnergy += trackster->regressed_energy();
0490 }
0491 cand.setRawEnergy(rawEnergy);
0492
0493 auto pdgId = cp.pdgId();
0494 auto charge = cp.charge();
0495 if (cand.trackPtr().isNonnull()) {
0496 auto const& track = cand.trackPtr().get();
0497 if (std::abs(pdgId) == 13) {
0498 cand.setPdgId(pdgId);
0499 } else {
0500 cand.setPdgId((isHad(pdgId) ? 211 : 11) * charge);
0501 }
0502 cand.setCharge(charge);
0503 math::XYZTLorentzVector p4(regressedEnergy * track->momentum().unit().x(),
0504 regressedEnergy * track->momentum().unit().y(),
0505 regressedEnergy * track->momentum().unit().z(),
0506 regressedEnergy);
0507 cand.setP4(p4);
0508 } else {
0509
0510
0511 if (charge != 0)
0512 cand.setPdgId(isHad(pdgId) ? 211 : 11);
0513 else if (pdgId == 111)
0514 cand.setPdgId(pdgId);
0515 else
0516 cand.setPdgId(isHad(pdgId) ? 130 : 22);
0517 cand.setCharge(0);
0518
0519 auto particleType = tracksterParticleTypeFromPdgId(cand.pdgId(), 1);
0520 cand.setIdProbability(particleType, 1.f);
0521
0522 const auto& simTracksterFromCP = (*result_fromCP)[i];
0523 float regressedEnergy = simTracksterFromCP.regressed_energy();
0524 math::XYZTLorentzVector p4(regressedEnergy * simTracksterFromCP.barycenter().unit().x(),
0525 regressedEnergy * simTracksterFromCP.barycenter().unit().y(),
0526 regressedEnergy * simTracksterFromCP.barycenter().unit().z(),
0527 regressedEnergy);
0528 cand.setP4(p4);
0529 }
0530 }
0531
0532 std::vector<int> all_nums(result_fromCP->size());
0533 std::iota(all_nums.begin(), all_nums.end(), 0);
0534
0535 std::vector<int> toRemove;
0536 std::set_difference(all_nums.begin(), all_nums.end(), toKeep.begin(), toKeep.end(), std::back_inserter(toRemove));
0537 std::sort(toRemove.begin(), toRemove.end(), [](int x, int y) { return x > y; });
0538 for (auto const& r : toRemove) {
0539 result_fromCP->erase(result_fromCP->begin() + r);
0540 result_ticlCandidates->erase(result_ticlCandidates->begin() + r);
0541 }
0542 evt.put(std::move(result_ticlCandidates));
0543 evt.put(std::move(output_mask));
0544 evt.put(std::move(result_fromCP), "fromCPs");
0545 evt.put(std::move(resultPU), "PU");
0546 evt.put(std::move(output_mask_fromCP), "fromCPs");
0547 evt.put(std::move(cpToSc_SimTrackstersMap));
0548 }