Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-06-22 02:23:59 

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/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       double 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);  // Conservative size, will call shrink_to_fit later
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     // Create a Trackster from the object entering HGCal
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     // Create a Trackster from any CP
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   // TODO: remove time computation from PCA calculation and
0367   //       store time from boundary position in simTracksters
0368   ticl::assignPCAtoTracksters(*result,
0369                               layerClusters,
0370                               layerClustersTimes,
0371                               rhtools_.getPositionLayer(rhtools_.lastLayerEE(doNose_)).z(),
0372                               computeLocalTime_);
0373   result->shrink_to_fit();
0374   ticl::assignPCAtoTracksters(*result_fromCP,
0375                               layerClusters,
0376                               layerClustersTimes,
0377                               rhtools_.getPositionLayer(rhtools_.lastLayerEE(doNose_)).z(),
0378                               computeLocalTime_);
0379 
0380   makePUTrackster(inputClusterMask, *output_mask, *resultPU, caloParticles_h.id(), 0);
0381 
0382   auto simTrackToRecoTrack = [&](UniqueSimTrackId simTkId) -> std::pair<int, float> {
0383     int trackIdx = -1;
0384     float quality = 0.f;
0385     auto ipos = simTrackToTPMap.mapping.find(simTkId);
0386     if (ipos != simTrackToTPMap.mapping.end()) {
0387       auto jpos = TPtoRecoTrackMap.find((ipos->second));
0388       if (jpos != TPtoRecoTrackMap.end()) {
0389         auto& associatedRecoTracks = jpos->val;
0390         if (!associatedRecoTracks.empty()) {
0391           // associated reco tracks are sorted by decreasing quality
0392           if (associatedRecoTracks[0].second > qualityCutTrack_) {
0393             trackIdx = &(*associatedRecoTracks[0].first) - &recoTracks[0];
0394             quality = associatedRecoTracks[0].second;
0395           }
0396         }
0397       }
0398     }
0399     return {trackIdx, quality};
0400   };
0401 
0402   // Creating the map from TrackingParticle to SimTrackstersFromCP
0403   auto& simTrackstersFromCP = *result_fromCP;
0404   for (unsigned int i = 0; i < simTrackstersFromCP.size(); ++i) {
0405     if (simTrackstersFromCP[i].vertices().empty())
0406       continue;
0407     const auto& simTrack = caloparticles[simTrackstersFromCP[i].seedIndex()].g4Tracks()[0];
0408     UniqueSimTrackId simTkIds(simTrack.trackId(), simTrack.eventId());
0409     auto bestAssociatedRecoTrack = simTrackToRecoTrack(simTkIds);
0410     if (bestAssociatedRecoTrack.first != -1 and bestAssociatedRecoTrack.second > qualityCutTrack_) {
0411       auto trackIndex = bestAssociatedRecoTrack.first;
0412       simTrackstersFromCP[i].setTrackIdx(trackIndex);
0413     }
0414   }
0415 
0416   auto& simTracksters = *result;
0417   // Creating the map from TrackingParticle to SimTrackster
0418   std::unordered_map<unsigned int, std::vector<unsigned int>> TPtoSimTracksterMap;
0419   for (unsigned int i = 0; i < simTracksters.size(); ++i) {
0420     const auto& simTrack = (simTracksters[i].seedID() == caloParticles_h.id())
0421                                ? caloparticles[simTracksters[i].seedIndex()].g4Tracks()[0]
0422                                : simclusters[simTracksters[i].seedIndex()].g4Tracks()[0];
0423     UniqueSimTrackId simTkIds(simTrack.trackId(), simTrack.eventId());
0424     auto bestAssociatedRecoTrack = simTrackToRecoTrack(simTkIds);
0425     if (bestAssociatedRecoTrack.first != -1 and bestAssociatedRecoTrack.second > qualityCutTrack_) {
0426       auto trackIndex = bestAssociatedRecoTrack.first;
0427       simTracksters[i].setTrackIdx(trackIndex);
0428     }
0429   }
0430 
0431   edm::OrphanHandle<std::vector<Trackster>> simTracksters_h = evt.put(std::move(result));
0432 
0433   // map between simTrack and Mtd SimTracksters to loop on them only one
0434   std::unordered_map<unsigned int, const MtdSimTrackster*> SimTrackToMtdST;
0435   for (unsigned int i = 0; i < MTDSimTracksters_h->size(); ++i) {
0436     const auto& simTrack = (*MTDSimTracksters_h)[i].g4Tracks()[0];
0437     SimTrackToMtdST[simTrack.trackId()] = &((*MTDSimTracksters_h)[i]);
0438   }
0439 
0440   result_ticlCandidates->resize(result_fromCP->size());
0441   std::vector<int> toKeep;
0442   for (size_t i = 0; i < simTracksters_h->size(); ++i) {
0443     const auto& simTrackster = (*simTracksters_h)[i];
0444     int cp_index = (simTrackster.seedID() == caloParticles_h.id())
0445                        ? simTrackster.seedIndex()
0446                        : SimClusterToCaloParticleMap[simTrackster.seedIndex()];
0447     auto const& tCP = (*result_fromCP)[cp_index];
0448     if (!tCP.vertices().empty()) {
0449       auto trackIndex = tCP.trackIdx();
0450 
0451       auto& cand = (*result_ticlCandidates)[cp_index];
0452       cand.addTrackster(edm::Ptr<Trackster>(simTracksters_h, i));
0453       if (trackIndex != -1 && caloparticles[cp_index].charge() != 0)
0454         cand.setTrackPtr(edm::Ptr<reco::Track>(recoTracks_h, trackIndex));
0455       toKeep.push_back(cp_index);
0456     }
0457   }
0458 
0459   auto isHad = [](int pdgId) {
0460     pdgId = std::abs(pdgId);
0461     if (pdgId == 111)
0462       return false;
0463     return (pdgId > 100 and pdgId < 900) or (pdgId > 1000 and pdgId < 9000);
0464   };
0465 
0466   for (size_t i = 0; i < result_ticlCandidates->size(); ++i) {
0467     auto cp_index = (*result_fromCP)[i].seedIndex();
0468     if (cp_index < 0)
0469       continue;
0470     auto& cand = (*result_ticlCandidates)[i];
0471     const auto& cp = caloparticles[cp_index];
0472     float rawEnergy = 0.f;
0473     float regressedEnergy = 0.f;
0474 
0475     const auto& simTrack = cp.g4Tracks()[0];
0476     auto pos = SimTrackToMtdST.find(simTrack.trackId());
0477     if (pos != SimTrackToMtdST.end()) {
0478       auto MTDst = pos->second;
0479       // TODO: once the associators have been implemented check if the MTDst is associated with a reco before adding the MTD time
0480       cand.setMTDTime(MTDst->time(), 0);
0481     }
0482 
0483     cand.setTime(simVertices[cp.g4Tracks()[0].vertIndex()].position().t() * pow(10, 9), 0);
0484 
0485     for (const auto& trackster : cand.tracksters()) {
0486       rawEnergy += trackster->raw_energy();
0487       regressedEnergy += trackster->regressed_energy();
0488     }
0489     cand.setRawEnergy(rawEnergy);
0490 
0491     auto pdgId = cp.pdgId();
0492     auto charge = cp.charge();
0493     if (cand.trackPtr().isNonnull()) {
0494       auto const& track = cand.trackPtr().get();
0495       if (std::abs(pdgId) == 13) {
0496         cand.setPdgId(pdgId);
0497       } else {
0498         cand.setPdgId((isHad(pdgId) ? 211 : 11) * charge);
0499       }
0500       cand.setCharge(charge);
0501       math::XYZTLorentzVector p4(regressedEnergy * track->momentum().unit().x(),
0502                                  regressedEnergy * track->momentum().unit().y(),
0503                                  regressedEnergy * track->momentum().unit().z(),
0504                                  regressedEnergy);
0505       cand.setP4(p4);
0506     } else {  // neutral candidates
0507       // a neutral candidate with a charged CaloParticle is charged without a reco track associated with it
0508       // set the charge = 0, but keep the real pdgId to keep track of that
0509       if (charge != 0)
0510         cand.setPdgId(isHad(pdgId) ? 211 : 11);
0511       else if (pdgId == 111)
0512         cand.setPdgId(pdgId);
0513       else
0514         cand.setPdgId(isHad(pdgId) ? 130 : 22);
0515       cand.setCharge(0);
0516 
0517       auto particleType = tracksterParticleTypeFromPdgId(cand.pdgId(), 1);
0518       cand.setIdProbability(particleType, 1.f);
0519 
0520       const auto& simTracksterFromCP = (*result_fromCP)[i];
0521       float regressedEnergy = simTracksterFromCP.regressed_energy();
0522       math::XYZTLorentzVector p4(regressedEnergy * simTracksterFromCP.barycenter().unit().x(),
0523                                  regressedEnergy * simTracksterFromCP.barycenter().unit().y(),
0524                                  regressedEnergy * simTracksterFromCP.barycenter().unit().z(),
0525                                  regressedEnergy);
0526       cand.setP4(p4);
0527     }
0528   }
0529 
0530   std::vector<int> all_nums(result_fromCP->size());  // vector containing all caloparticles indexes
0531   std::iota(all_nums.begin(), all_nums.end(), 0);    // fill the vector with consecutive numbers starting from 0
0532 
0533   std::vector<int> toRemove;
0534   std::set_difference(all_nums.begin(), all_nums.end(), toKeep.begin(), toKeep.end(), std::back_inserter(toRemove));
0535   std::sort(toRemove.begin(), toRemove.end(), [](int x, int y) { return x > y; });
0536   for (auto const& r : toRemove) {
0537     result_fromCP->erase(result_fromCP->begin() + r);
0538     result_ticlCandidates->erase(result_ticlCandidates->begin() + r);
0539   }
0540   evt.put(std::move(result_ticlCandidates));
0541   evt.put(std::move(output_mask));
0542   evt.put(std::move(result_fromCP), "fromCPs");
0543   evt.put(std::move(resultPU), "PU");
0544   evt.put(std::move(output_mask_fromCP), "fromCPs");
0545   evt.put(std::move(cpToSc_SimTrackstersMap));
0546 }
[ 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!