File indexing completed on 2024-09-11 04:33:30
0001
0002
0003
0004 #include "FWCore/Framework/interface/Frameworkfwd.h"
0005 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0006
0007 #include "FWCore/Framework/interface/ESHandle.h"
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/MakerMacros.h"
0010
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012
0013 #include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
0014 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0015 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0016 #include "DataFormats/HcalRecHit/interface/HFRecHit.h"
0017 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0018 #include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h"
0019 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0020 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0021 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0022 #include "DataFormats/VertexReco/interface/Vertex.h"
0023 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0024 #include "DataFormats/ParticleFlowReco/interface/PFBlockElementSuperCluster.h"
0025 #include "DataFormats/ParticleFlowReco/interface/PFBlockElementGsfTrack.h"
0026 #include "DataFormats/ParticleFlowReco/interface/PFBlockElementTrack.h"
0027 #include "DataFormats/ParticleFlowReco/interface/PFBlockElementBrem.h"
0028 #include "DataFormats/ParticleFlowReco/interface/PFBlockElementCluster.h"
0029 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0030 #include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
0031 #include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
0032 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0033 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0034 #include "SimDataFormats/Associations/interface/TrackAssociation.h"
0035 #include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
0036 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0037 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
0038
0039 #include "DataFormats/Math/interface/deltaPhi.h"
0040 #include "DataFormats/TrackReco/interface/Track.h"
0041 #include "DataFormats/DetId/interface/DetId.h"
0042
0043 #include "DataFormats/GeometrySurface/interface/PlaneBuilder.h"
0044 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0045 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0046 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0047 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0048
0049 #include "CommonTools/BaseParticlePropagator/interface/BaseParticlePropagator.h"
0050 #include "Math/Transform3D.h"
0051 #include "MagneticField/Engine/interface/MagneticField.h"
0052 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0053 #include "MagneticField/VolumeGeometry/interface/MagVolumeOutsideValidity.h"
0054 #include "RecoParticleFlow/PFProducer/interface/MLPFModel.h"
0055
0056 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0057 #include "FWCore/ServiceRegistry/interface/Service.h"
0058 #include "TH1F.h"
0059 #include "TVector2.h"
0060 #include "TTree.h"
0061
0062 #include <map>
0063 #include <set>
0064 #include <string>
0065 #include <vector>
0066 #include <set>
0067
0068 using namespace std;
0069
0070 class ElementWithIndex {
0071 public:
0072 const reco::PFBlockElement& orig;
0073 size_t idx_block;
0074 size_t idx_elem;
0075 ElementWithIndex(const reco::PFBlockElement& _orig, size_t _idx_block, size_t _idx_elem)
0076 : orig(_orig), idx_block(_idx_block), idx_elem(_idx_elem) {}
0077 };
0078
0079 vector<int> find_element_ref(const vector<ElementWithIndex>& vec, const edm::RefToBase<reco::Track>& r) {
0080 vector<int> ret;
0081 for (unsigned int i = 0; i < vec.size(); i++) {
0082 const auto& elem = vec.at(i);
0083 if (elem.orig.type() == reco::PFBlockElement::TRACK) {
0084 const auto& ref = elem.orig.trackRef();
0085 if (ref.isNonnull() && ref->extra().isNonnull()) {
0086 if (ref.key() == r.key()) {
0087 ret.push_back(i);
0088 }
0089 }
0090 } else if (elem.orig.type() == reco::PFBlockElement::GSF) {
0091 const auto& ref = ((const reco::PFBlockElementGsfTrack*)&elem.orig)->clusterRef();
0092 if (ref.isNonnull()) {
0093 if (ref.key() == r.key()) {
0094 ret.push_back(i);
0095 }
0096 }
0097 } else if (elem.orig.type() == reco::PFBlockElement::BREM) {
0098 const auto& ref = ((const reco::PFBlockElementBrem*)&elem.orig)->clusterRef();
0099 if (ref.isNonnull()) {
0100 if (ref.key() == r.key()) {
0101 ret.push_back(i);
0102 }
0103 }
0104 }
0105 }
0106 return ret;
0107 }
0108
0109 double detid_compare(const map<uint64_t, double>& rechits, const map<uint64_t, double>& simhits) {
0110 double ret = 0.0;
0111
0112 for (const auto& rh : rechits) {
0113 for (const auto& sh : simhits) {
0114 if (rh.first == sh.first) {
0115
0116 ret += rh.second * sh.second;
0117 break;
0118 }
0119 }
0120 }
0121 return ret;
0122 }
0123
0124 class PFAnalysis : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0125 public:
0126 typedef ROOT::Math::Transform3D::Point Point;
0127
0128 PFAnalysis();
0129 explicit PFAnalysis(const edm::ParameterSet&);
0130 ~PFAnalysis() override;
0131
0132 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0133 void beginRun(edm::Run const& iEvent, edm::EventSetup const&) override;
0134 void endRun(edm::Run const& iEvent, edm::EventSetup const&) override;
0135
0136 private:
0137 void beginJob() override;
0138 void analyze(const edm::Event&, const edm::EventSetup&) override;
0139 void endJob() override;
0140
0141 void processTrackingParticles(const edm::View<TrackingParticle>& trackingParticles,
0142 edm::Handle<edm::View<TrackingParticle>>& trackingParticlesHandle);
0143
0144 pair<vector<ElementWithIndex>, vector<tuple<int, int, float>>> processBlocks(
0145 const std::vector<reco::PFBlock>& pfBlocks);
0146
0147 void associateClusterToSimCluster(const vector<ElementWithIndex>& all_elements);
0148
0149 void clearVariables();
0150
0151 GlobalPoint getHitPosition(const DetId& id);
0152
0153
0154 edm::EDGetTokenT<std::vector<reco::GenParticle>> genParticles_;
0155 edm::EDGetTokenT<edm::View<TrackingParticle>> trackingParticles_;
0156 edm::EDGetTokenT<edm::View<CaloParticle>> caloParticles_;
0157 edm::EDGetTokenT<edm::View<reco::Track>> tracks_;
0158 edm::EDGetTokenT<std::vector<reco::PFBlock>> pfBlocks_;
0159 edm::EDGetTokenT<std::vector<reco::PFCandidate>> pfCandidates_;
0160 edm::EDGetTokenT<reco::RecoToSimCollection> tracks_recotosim_;
0161
0162 TTree* t_;
0163
0164 edm::RunNumber_t ev_run_;
0165 edm::LuminosityBlockNumber_t ev_lumi_;
0166 edm::EventNumber_t ev_event_;
0167
0168 vector<float> trackingparticle_eta_;
0169 vector<float> trackingparticle_phi_;
0170 vector<float> trackingparticle_pt_;
0171 vector<float> trackingparticle_px_;
0172 vector<float> trackingparticle_py_;
0173 vector<float> trackingparticle_pz_;
0174 vector<float> trackingparticle_energy_;
0175 vector<float> trackingparticle_dvx_;
0176 vector<float> trackingparticle_dvy_;
0177 vector<float> trackingparticle_dvz_;
0178 vector<int> trackingparticle_bx_;
0179 vector<int> trackingparticle_ev_;
0180 vector<float> trackingparticle_ovx_;
0181 vector<float> trackingparticle_ovy_;
0182 vector<float> trackingparticle_ovz_;
0183 vector<float> trackingparticle_exx_;
0184 vector<float> trackingparticle_exy_;
0185 vector<int> trackingparticle_mother_;
0186 vector<int> trackingparticle_pid_;
0187
0188 vector<float> simcluster_eta_;
0189 vector<float> simcluster_phi_;
0190 vector<float> simcluster_pt_;
0191 vector<float> simcluster_energy_;
0192 vector<float> simcluster_px_;
0193 vector<float> simcluster_py_;
0194 vector<float> simcluster_pz_;
0195 vector<int> simcluster_bx_;
0196 vector<int> simcluster_ev_;
0197 vector<int> simcluster_pid_;
0198 vector<int> simcluster_idx_trackingparticle_;
0199 vector<int> simcluster_nhits_;
0200 vector<std::map<uint64_t, double>> simcluster_detids_;
0201
0202 vector<float> simhit_frac_;
0203 vector<float> simhit_x_;
0204 vector<float> simhit_y_;
0205 vector<float> simhit_z_;
0206 vector<float> simhit_eta_;
0207 vector<float> simhit_phi_;
0208 vector<int> simhit_det_;
0209 vector<int> simhit_subdet_;
0210 vector<int> simhit_idx_simcluster_;
0211 vector<uint64_t> simhit_detid_;
0212
0213 vector<float> rechit_e_;
0214 vector<float> rechit_x_;
0215 vector<float> rechit_y_;
0216 vector<float> rechit_z_;
0217 vector<float> rechit_det_;
0218 vector<float> rechit_subdet_;
0219 vector<float> rechit_eta_;
0220 vector<float> rechit_phi_;
0221 vector<int> rechit_idx_element_;
0222 vector<uint64_t> rechit_detid_;
0223
0224 vector<float> simtrack_x_;
0225 vector<float> simtrack_y_;
0226 vector<float> simtrack_z_;
0227 vector<int> simtrack_idx_simcluster_;
0228 vector<int> simtrack_pid_;
0229
0230 vector<float> gen_eta_;
0231 vector<float> gen_phi_;
0232 vector<float> gen_pt_;
0233 vector<float> gen_px_;
0234 vector<float> gen_py_;
0235 vector<float> gen_pz_;
0236 vector<float> gen_energy_;
0237 vector<int> gen_charge_;
0238 vector<int> gen_pdgid_;
0239 vector<int> gen_status_;
0240 vector<vector<int>> gen_daughters_;
0241
0242 vector<float> element_pt_;
0243 vector<float> element_px_;
0244 vector<float> element_py_;
0245 vector<float> element_pz_;
0246 vector<float> element_deltap_;
0247 vector<float> element_sigmadeltap_;
0248 vector<float> element_eta_;
0249 vector<float> element_phi_;
0250 vector<float> element_energy_;
0251 vector<float> element_corr_energy_;
0252 vector<float> element_eta_ecal_;
0253 vector<float> element_phi_ecal_;
0254 vector<float> element_eta_hcal_;
0255 vector<float> element_phi_hcal_;
0256 vector<int> element_charge_;
0257 vector<int> element_type_;
0258 vector<int> element_layer_;
0259 vector<float> element_depth_;
0260 vector<float> element_trajpoint_;
0261 vector<float> element_muon_dt_hits_;
0262 vector<float> element_muon_csc_hits_;
0263 vector<float> element_muon_type_;
0264 vector<float> element_cluster_flags_;
0265 vector<float> element_gsf_electronseed_trkorecal_;
0266 vector<float> element_num_hits_;
0267
0268 vector<int> element_distance_i_;
0269 vector<int> element_distance_j_;
0270 vector<float> element_distance_d_;
0271
0272 vector<float> pfcandidate_eta_;
0273 vector<float> pfcandidate_phi_;
0274 vector<float> pfcandidate_pt_;
0275 vector<float> pfcandidate_px_;
0276 vector<float> pfcandidate_py_;
0277 vector<float> pfcandidate_pz_;
0278 vector<float> pfcandidate_energy_;
0279 vector<int> pfcandidate_pdgid_;
0280
0281 vector<pair<int, int>> trackingparticle_to_element;
0282 vector<pair<int, int>> simcluster_to_element;
0283 vector<float> simcluster_to_element_cmp;
0284 vector<pair<int, int>> element_to_candidate;
0285
0286
0287 MagneticField const* aField_;
0288
0289 edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geometryToken_;
0290 edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> topologyToken_;
0291 edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken_;
0292 edm::ESGetToken<HcalDDDRecConstants, HcalRecNumberingRecord> hcalDDDrecToken_;
0293
0294 CaloGeometry* geom;
0295 HcalTopology* hcal_topo;
0296 const HcalDDDRecConstants* hcons;
0297
0298 bool saveHits;
0299 };
0300
0301 PFAnalysis::PFAnalysis() { ; }
0302
0303 PFAnalysis::PFAnalysis(const edm::ParameterSet& iConfig) {
0304 tracks_recotosim_ = consumes<reco::RecoToSimCollection>(edm::InputTag("trackingParticleRecoTrackAsssociation"));
0305 trackingParticles_ = consumes<edm::View<TrackingParticle>>(edm::InputTag("mix", "MergedTrackTruth"));
0306 caloParticles_ = consumes<edm::View<CaloParticle>>(edm::InputTag("mix", "MergedCaloTruth"));
0307 genParticles_ = consumes<std::vector<reco::GenParticle>>(edm::InputTag("genParticles"));
0308 pfBlocks_ = consumes<std::vector<reco::PFBlock>>(edm::InputTag("particleFlowBlock"));
0309 pfCandidates_ = consumes<std::vector<reco::PFCandidate>>(edm::InputTag("particleFlow"));
0310 tracks_ = consumes<edm::View<reco::Track>>(edm::InputTag("generalTracks"));
0311 saveHits = iConfig.getUntrackedParameter<bool>("saveHits", false);
0312
0313 geometryToken_ = esConsumes<CaloGeometry, CaloGeometryRecord>(edm::ESInputTag{});
0314 topologyToken_ = esConsumes<HcalTopology, HcalRecNumberingRecord>(edm::ESInputTag{});
0315 magFieldToken_ = esConsumes<edm::Transition::BeginRun>();
0316 hcalDDDrecToken_ = esConsumes<edm::Transition::BeginRun>();
0317
0318 usesResource(TFileService::kSharedResource);
0319 edm::Service<TFileService> fs;
0320 fs->make<TH1F>("total", "total", 100, 0, 5.);
0321
0322 t_ = fs->make<TTree>("pftree", "pftree");
0323
0324
0325 t_->Branch("event", &ev_event_);
0326 t_->Branch("lumi", &ev_lumi_);
0327 t_->Branch("run", &ev_run_);
0328
0329 t_->Branch("trackingparticle_eta", &trackingparticle_eta_);
0330 t_->Branch("trackingparticle_phi", &trackingparticle_phi_);
0331 t_->Branch("trackingparticle_pt", &trackingparticle_pt_);
0332 t_->Branch("trackingparticle_px", &trackingparticle_px_);
0333 t_->Branch("trackingparticle_py", &trackingparticle_py_);
0334 t_->Branch("trackingparticle_pz", &trackingparticle_pz_);
0335 t_->Branch("trackingparticle_energy", &trackingparticle_energy_);
0336 t_->Branch("trackingparticle_dvx", &trackingparticle_dvx_);
0337 t_->Branch("trackingparticle_dvy", &trackingparticle_dvy_);
0338 t_->Branch("trackingparticle_dvz", &trackingparticle_dvz_);
0339 t_->Branch("trackingparticle_bx", &trackingparticle_bx_);
0340 t_->Branch("trackingparticle_ev", &trackingparticle_ev_);
0341 t_->Branch("trackingparticle_pid", &trackingparticle_pid_);
0342
0343 t_->Branch("simcluster_eta", &simcluster_eta_);
0344 t_->Branch("simcluster_phi", &simcluster_phi_);
0345 t_->Branch("simcluster_pt", &simcluster_pt_);
0346 t_->Branch("simcluster_px", &simcluster_px_);
0347 t_->Branch("simcluster_py", &simcluster_py_);
0348 t_->Branch("simcluster_pz", &simcluster_pz_);
0349 t_->Branch("simcluster_energy", &simcluster_energy_);
0350 t_->Branch("simcluster_bx", &simcluster_bx_);
0351 t_->Branch("simcluster_ev", &simcluster_ev_);
0352 t_->Branch("simcluster_pid", &simcluster_pid_);
0353 t_->Branch("simcluster_idx_trackingparticle", &simcluster_idx_trackingparticle_);
0354 t_->Branch("simcluster_nhits", &simcluster_nhits_);
0355
0356 if (saveHits) {
0357 t_->Branch("simhit_frac", &simhit_frac_);
0358 t_->Branch("simhit_x", &simhit_x_);
0359 t_->Branch("simhit_y", &simhit_y_);
0360 t_->Branch("simhit_z", &simhit_z_);
0361 t_->Branch("simhit_det", &simhit_det_);
0362 t_->Branch("simhit_subdet", &simhit_subdet_);
0363 t_->Branch("simhit_eta", &simhit_eta_);
0364 t_->Branch("simhit_phi", &simhit_phi_);
0365 t_->Branch("simhit_idx_simcluster", &simhit_idx_simcluster_);
0366 t_->Branch("simhit_detid", &simhit_detid_);
0367
0368 t_->Branch("rechit_e", &rechit_e_);
0369 t_->Branch("rechit_x", &rechit_x_);
0370 t_->Branch("rechit_y", &rechit_y_);
0371 t_->Branch("rechit_z", &rechit_z_);
0372 t_->Branch("rechit_det", &rechit_det_);
0373 t_->Branch("rechit_subdet", &rechit_subdet_);
0374 t_->Branch("rechit_eta", &rechit_eta_);
0375 t_->Branch("rechit_phi", &rechit_phi_);
0376 t_->Branch("rechit_idx_element", &rechit_idx_element_);
0377 t_->Branch("rechit_detid", &rechit_detid_);
0378 }
0379
0380 t_->Branch("simtrack_x", &simtrack_x_);
0381 t_->Branch("simtrack_y", &simtrack_y_);
0382 t_->Branch("simtrack_z", &simtrack_z_);
0383 t_->Branch("simtrack_idx_simcluster_", &simtrack_idx_simcluster_);
0384 t_->Branch("simtrack_pid", &simtrack_pid_);
0385
0386 t_->Branch("gen_eta", &gen_eta_);
0387 t_->Branch("gen_phi", &gen_phi_);
0388 t_->Branch("gen_pt", &gen_pt_);
0389 t_->Branch("gen_px", &gen_px_);
0390 t_->Branch("gen_py", &gen_py_);
0391 t_->Branch("gen_pz", &gen_pz_);
0392 t_->Branch("gen_energy", &gen_energy_);
0393 t_->Branch("gen_charge", &gen_charge_);
0394 t_->Branch("gen_pdgid", &gen_pdgid_);
0395 t_->Branch("gen_status", &gen_status_);
0396 t_->Branch("gen_daughters", &gen_daughters_);
0397
0398
0399 t_->Branch("element_pt", &element_pt_);
0400 t_->Branch("element_px", &element_px_);
0401 t_->Branch("element_py", &element_py_);
0402 t_->Branch("element_pz", &element_pz_);
0403 t_->Branch("element_deltap", &element_deltap_);
0404 t_->Branch("element_sigmadeltap", &element_sigmadeltap_);
0405 t_->Branch("element_eta", &element_eta_);
0406 t_->Branch("element_phi", &element_phi_);
0407 t_->Branch("element_energy", &element_energy_);
0408 t_->Branch("element_corr_energy", &element_corr_energy_);
0409 t_->Branch("element_eta_ecal", &element_eta_ecal_);
0410 t_->Branch("element_phi_ecal", &element_phi_ecal_);
0411 t_->Branch("element_eta_hcal", &element_eta_hcal_);
0412 t_->Branch("element_phi_hcal", &element_phi_hcal_);
0413 t_->Branch("element_charge", &element_charge_);
0414 t_->Branch("element_type", &element_type_);
0415 t_->Branch("element_layer", &element_layer_);
0416 t_->Branch("element_depth", &element_depth_);
0417 t_->Branch("element_trajpoint", &element_trajpoint_);
0418 t_->Branch("element_muon_dt_hits", &element_muon_dt_hits_);
0419 t_->Branch("element_muon_csc_hits", &element_muon_csc_hits_);
0420 t_->Branch("element_muon_type", &element_muon_type_);
0421 t_->Branch("element_cluster_flags", &element_cluster_flags_);
0422 t_->Branch("element_gsf_electronseed_trkorecal", &element_gsf_electronseed_trkorecal_);
0423 t_->Branch("element_num_hits", &element_num_hits_);
0424
0425
0426 t_->Branch("element_distance_i", &element_distance_i_);
0427 t_->Branch("element_distance_j", &element_distance_j_);
0428 t_->Branch("element_distance_d", &element_distance_d_);
0429
0430 t_->Branch("pfcandidate_eta", &pfcandidate_eta_);
0431 t_->Branch("pfcandidate_phi", &pfcandidate_phi_);
0432 t_->Branch("pfcandidate_pt", &pfcandidate_pt_);
0433 t_->Branch("pfcandidate_px", &pfcandidate_px_);
0434 t_->Branch("pfcandidate_py", &pfcandidate_py_);
0435 t_->Branch("pfcandidate_pz", &pfcandidate_pz_);
0436 t_->Branch("pfcandidate_energy", &pfcandidate_energy_);
0437 t_->Branch("pfcandidate_pdgid", &pfcandidate_pdgid_);
0438
0439
0440 t_->Branch("trackingparticle_to_element", &trackingparticle_to_element);
0441 t_->Branch("simcluster_to_element", &simcluster_to_element);
0442 t_->Branch("simcluster_to_element_cmp", &simcluster_to_element_cmp);
0443 t_->Branch("element_to_candidate", &element_to_candidate);
0444 }
0445
0446 PFAnalysis::~PFAnalysis() {}
0447
0448 void PFAnalysis::clearVariables() {
0449 ev_run_ = 0;
0450 ev_lumi_ = 0;
0451 ev_event_ = 0;
0452
0453 trackingparticle_to_element.clear();
0454 simcluster_to_element.clear();
0455 simcluster_to_element_cmp.clear();
0456 element_to_candidate.clear();
0457
0458 trackingparticle_eta_.clear();
0459 trackingparticle_phi_.clear();
0460 trackingparticle_pt_.clear();
0461 trackingparticle_px_.clear();
0462 trackingparticle_py_.clear();
0463 trackingparticle_pz_.clear();
0464 trackingparticle_energy_.clear();
0465 trackingparticle_dvx_.clear();
0466 trackingparticle_dvy_.clear();
0467 trackingparticle_dvz_.clear();
0468 trackingparticle_bx_.clear();
0469 trackingparticle_ev_.clear();
0470 trackingparticle_ovx_.clear();
0471 trackingparticle_ovy_.clear();
0472 trackingparticle_ovz_.clear();
0473 trackingparticle_exx_.clear();
0474 trackingparticle_exy_.clear();
0475 trackingparticle_mother_.clear();
0476 trackingparticle_pid_.clear();
0477
0478 simcluster_eta_.clear();
0479 simcluster_phi_.clear();
0480 simcluster_pt_.clear();
0481 simcluster_energy_.clear();
0482 simcluster_pid_.clear();
0483 simcluster_detids_.clear();
0484 simcluster_bx_.clear();
0485 simcluster_ev_.clear();
0486 simcluster_px_.clear();
0487 simcluster_py_.clear();
0488 simcluster_pz_.clear();
0489 simcluster_idx_trackingparticle_.clear();
0490 simcluster_nhits_.clear();
0491
0492 if (saveHits) {
0493 simhit_frac_.clear();
0494 simhit_x_.clear();
0495 simhit_y_.clear();
0496 simhit_z_.clear();
0497 simhit_det_.clear();
0498 simhit_subdet_.clear();
0499 simhit_eta_.clear();
0500 simhit_phi_.clear();
0501 simhit_idx_simcluster_.clear();
0502 simhit_detid_.clear();
0503
0504 rechit_e_.clear();
0505 rechit_x_.clear();
0506 rechit_y_.clear();
0507 rechit_z_.clear();
0508 rechit_det_.clear();
0509 rechit_subdet_.clear();
0510 rechit_eta_.clear();
0511 rechit_phi_.clear();
0512 rechit_idx_element_.clear();
0513 rechit_detid_.clear();
0514 }
0515
0516 simtrack_x_.clear();
0517 simtrack_y_.clear();
0518 simtrack_z_.clear();
0519 simtrack_idx_simcluster_.clear();
0520 simtrack_pid_.clear();
0521
0522 gen_eta_.clear();
0523 gen_phi_.clear();
0524 gen_pt_.clear();
0525 gen_px_.clear();
0526 gen_py_.clear();
0527 gen_pz_.clear();
0528 gen_energy_.clear();
0529 gen_charge_.clear();
0530 gen_pdgid_.clear();
0531 gen_status_.clear();
0532 gen_daughters_.clear();
0533
0534 element_pt_.clear();
0535 element_px_.clear();
0536 element_py_.clear();
0537 element_pz_.clear();
0538 element_deltap_.clear();
0539 element_sigmadeltap_.clear();
0540 element_eta_.clear();
0541 element_phi_.clear();
0542 element_energy_.clear();
0543 element_corr_energy_.clear();
0544 element_eta_ecal_.clear();
0545 element_phi_ecal_.clear();
0546 element_eta_hcal_.clear();
0547 element_phi_hcal_.clear();
0548 element_charge_.clear();
0549 element_type_.clear();
0550 element_layer_.clear();
0551 element_depth_.clear();
0552 element_trajpoint_.clear();
0553 element_muon_dt_hits_.clear();
0554 element_muon_csc_hits_.clear();
0555 element_muon_type_.clear();
0556 element_cluster_flags_.clear();
0557 element_gsf_electronseed_trkorecal_.clear();
0558 element_num_hits_.clear();
0559
0560 element_distance_i_.clear();
0561 element_distance_j_.clear();
0562 element_distance_d_.clear();
0563
0564 pfcandidate_eta_.clear();
0565 pfcandidate_phi_.clear();
0566 pfcandidate_pt_.clear();
0567 pfcandidate_px_.clear();
0568 pfcandidate_py_.clear();
0569 pfcandidate_pz_.clear();
0570 pfcandidate_energy_.clear();
0571 pfcandidate_pdgid_.clear();
0572
0573 }
0574
0575 GlobalPoint PFAnalysis::getHitPosition(const DetId& id) {
0576 GlobalPoint ret;
0577
0578 bool present = false;
0579 if (((id.det() == DetId::Ecal &&
0580 (id.subdetId() == EcalBarrel || id.subdetId() == EcalEndcap || id.subdetId() == EcalPreshower)) ||
0581 (id.det() == DetId::Hcal && (id.subdetId() == HcalBarrel || id.subdetId() == HcalEndcap ||
0582 id.subdetId() == HcalForward || id.subdetId() == HcalOuter)))) {
0583 const CaloSubdetectorGeometry* geom_sd(geom->getSubdetectorGeometry(id.det(), id.subdetId()));
0584 present = geom_sd->present(id);
0585 if (present) {
0586 const auto& cell = geom_sd->getGeometry(id);
0587 ret = GlobalPoint(cell->getPosition());
0588 }
0589 }
0590 return ret;
0591 }
0592
0593 void PFAnalysis::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0594 clearVariables();
0595
0596 auto& pG = iSetup.getData(geometryToken_);
0597 geom = (CaloGeometry*)&pG;
0598 auto& pT = iSetup.getData(topologyToken_);
0599 hcal_topo = (HcalTopology*)&pT;
0600
0601
0602 edm::Handle<edm::View<TrackingParticle>> trackingParticlesHandle;
0603 iEvent.getByToken(trackingParticles_, trackingParticlesHandle);
0604 const edm::View<TrackingParticle>& trackingParticles = *trackingParticlesHandle;
0605
0606 edm::Handle<edm::View<CaloParticle>> caloParticlesHandle;
0607 iEvent.getByToken(caloParticles_, caloParticlesHandle);
0608 const edm::View<CaloParticle>& caloParticles = *caloParticlesHandle;
0609
0610
0611 edm::Handle<reco::RecoToSimCollection> recotosimCollection;
0612 iEvent.getByToken(tracks_recotosim_, recotosimCollection);
0613 const auto recotosim = *recotosimCollection;
0614
0615 edm::Handle<edm::View<reco::Track>> trackHandle;
0616 iEvent.getByToken(tracks_, trackHandle);
0617 const edm::View<reco::Track>& tracks = *trackHandle;
0618
0619 edm::Handle<std::vector<reco::GenParticle>> genParticlesHandle;
0620 iEvent.getByToken(genParticles_, genParticlesHandle);
0621 for (std::vector<reco::GenParticle>::const_iterator it_p = genParticlesHandle->begin();
0622 it_p != genParticlesHandle->end();
0623 ++it_p) {
0624 gen_eta_.push_back(it_p->eta());
0625 gen_phi_.push_back(it_p->phi());
0626 gen_pt_.push_back(it_p->pt());
0627 gen_px_.push_back(it_p->px());
0628 gen_py_.push_back(it_p->py());
0629 gen_pz_.push_back(it_p->pz());
0630 gen_energy_.push_back(it_p->energy());
0631 gen_charge_.push_back(it_p->charge());
0632 gen_pdgid_.push_back(it_p->pdgId());
0633 gen_status_.push_back(it_p->status());
0634 std::vector<int> daughters(it_p->daughterRefVector().size(), 0);
0635 for (unsigned j = 0; j < it_p->daughterRefVector().size(); ++j) {
0636 daughters[j] = static_cast<int>(it_p->daughterRefVector().at(j).key());
0637 }
0638 gen_daughters_.push_back(daughters);
0639 }
0640
0641 edm::Handle<std::vector<reco::PFCandidate>> pfCandidatesHandle;
0642 iEvent.getByToken(pfCandidates_, pfCandidatesHandle);
0643 std::vector<reco::PFCandidate> pfCandidates = *pfCandidatesHandle;
0644
0645 edm::Handle<std::vector<reco::PFBlock>> pfBlocksHandle;
0646 iEvent.getByToken(pfBlocks_, pfBlocksHandle);
0647 std::vector<reco::PFBlock> pfBlocks = *pfBlocksHandle;
0648
0649
0650 const auto& all_elements_distances = processBlocks(pfBlocks);
0651 const auto& all_elements = all_elements_distances.first;
0652 const auto& all_distances = all_elements_distances.second;
0653 assert(!all_elements.empty());
0654
0655 for (const auto& d : all_distances) {
0656 element_distance_i_.push_back(get<0>(d));
0657 element_distance_j_.push_back(get<1>(d));
0658 element_distance_d_.push_back(get<2>(d));
0659 }
0660
0661
0662 for (unsigned long ntrack = 0; ntrack < tracks.size(); ntrack++) {
0663 edm::RefToBase<reco::Track> trackref(trackHandle, ntrack);
0664 const auto vec_idx_in_all_elements = find_element_ref(all_elements, trackref);
0665
0666
0667 if (vec_idx_in_all_elements.empty()) {
0668 continue;
0669 }
0670
0671 if (recotosim.find(trackref) != recotosim.end()) {
0672 const auto& tps = recotosim[trackref];
0673 for (const auto& tp : tps) {
0674 edm::Ref<std::vector<TrackingParticle>> tpr = tp.first;
0675 for (auto idx_in_all_elements : vec_idx_in_all_elements) {
0676 trackingparticle_to_element.emplace_back(tpr.key(), idx_in_all_elements);
0677 }
0678 }
0679 }
0680 }
0681
0682 processTrackingParticles(trackingParticles, trackingParticlesHandle);
0683
0684 int idx_simcluster = 0;
0685
0686 for (unsigned long ncaloparticle = 0; ncaloparticle < caloParticles.size(); ncaloparticle++) {
0687 const auto& cp = caloParticles.at(ncaloparticle);
0688 edm::RefToBase<CaloParticle> cpref(caloParticlesHandle, ncaloparticle);
0689
0690 int nhits = 0;
0691 for (const auto& simcluster : cp.simClusters()) {
0692
0693 map<uint64_t, double> detid_energy;
0694
0695 simcluster_nhits_.push_back(nhits);
0696 simcluster_eta_.push_back(simcluster->p4().eta());
0697 simcluster_phi_.push_back(simcluster->p4().phi());
0698 simcluster_pt_.push_back(simcluster->p4().pt());
0699 simcluster_energy_.push_back(simcluster->energy());
0700 simcluster_pid_.push_back(simcluster->pdgId());
0701 simcluster_bx_.push_back(simcluster->eventId().bunchCrossing());
0702 simcluster_ev_.push_back(simcluster->eventId().event());
0703
0704 simcluster_px_.push_back(simcluster->p4().x());
0705 simcluster_py_.push_back(simcluster->p4().y());
0706 simcluster_pz_.push_back(simcluster->p4().z());
0707
0708 for (const auto& hf : simcluster->hits_and_fractions()) {
0709 DetId id(hf.first);
0710
0711 if (id.det() == DetId::Hcal || id.det() == DetId::Ecal) {
0712 const auto& pos = getHitPosition(id);
0713 nhits += 1;
0714
0715 const float x = pos.x();
0716 const float y = pos.y();
0717 const float z = pos.z();
0718 const float eta = pos.eta();
0719 const float phi = pos.phi();
0720
0721 simhit_frac_.push_back(hf.second);
0722 simhit_x_.push_back(x);
0723 simhit_y_.push_back(y);
0724 simhit_z_.push_back(z);
0725 simhit_det_.push_back(id.det());
0726 simhit_subdet_.push_back(id.subdetId());
0727 simhit_eta_.push_back(eta);
0728 simhit_phi_.push_back(phi);
0729 simhit_idx_simcluster_.push_back(idx_simcluster);
0730 simhit_detid_.push_back(id.rawId());
0731 detid_energy[id.rawId()] += hf.second;
0732 }
0733 }
0734
0735 int simcluster_to_trackingparticle = -1;
0736 for (const auto& simtrack : simcluster->g4Tracks()) {
0737 simtrack_x_.push_back(simtrack.trackerSurfacePosition().x());
0738 simtrack_y_.push_back(simtrack.trackerSurfacePosition().y());
0739 simtrack_z_.push_back(simtrack.trackerSurfacePosition().z());
0740 simtrack_idx_simcluster_.push_back(idx_simcluster);
0741 simtrack_pid_.push_back(simtrack.type());
0742
0743 for (unsigned int itp = 0; itp < trackingParticles.size(); itp++) {
0744 const auto& simtrack2 = trackingParticles.at(itp).g4Tracks().at(0);
0745
0746 if (simtrack.eventId() == simtrack2.eventId() && simtrack.trackId() == simtrack2.trackId()) {
0747 simcluster_to_trackingparticle = itp;
0748
0749 break;
0750 }
0751 }
0752 }
0753
0754 simcluster_detids_.push_back(detid_energy);
0755 simcluster_idx_trackingparticle_.push_back(simcluster_to_trackingparticle);
0756
0757 idx_simcluster += 1;
0758 }
0759 }
0760
0761 associateClusterToSimCluster(all_elements);
0762
0763
0764 for (unsigned int ielem = 0; ielem < all_elements.size(); ielem++) {
0765 const auto& elem = all_elements.at(ielem);
0766 const auto& orig = elem.orig;
0767 reco::PFBlockElement::Type type = orig.type();
0768
0769 float pt = 0.0;
0770 float deltap = 0.0;
0771 float sigmadeltap = 0.0;
0772 float px = 0.0;
0773 float py = 0.0;
0774 float pz = 0.0;
0775 float eta = 0.0;
0776 float phi = 0.0;
0777 float energy = 0.0;
0778 float corr_energy = 0.0;
0779 float trajpoint = 0.0;
0780 float eta_ecal = 0.0;
0781 float phi_ecal = 0.0;
0782 float eta_hcal = 0.0;
0783 float phi_hcal = 0.0;
0784 int charge = 0;
0785 int layer = 0;
0786 float depth = 0;
0787 float muon_dt_hits = 0.0;
0788 float muon_csc_hits = 0.0;
0789 float muon_type = 0.0;
0790 float cluster_flags = 0.0;
0791 float gsf_electronseed_trkorecal = 0.0;
0792 float num_hits = 0.0;
0793
0794 if (type == reco::PFBlockElement::TRACK) {
0795 const auto& matched_pftrack = orig.trackRefPF();
0796 if (matched_pftrack.isNonnull()) {
0797 const auto& atECAL = matched_pftrack->extrapolatedPoint(reco::PFTrajectoryPoint::ECALShowerMax);
0798 const auto& atHCAL = matched_pftrack->extrapolatedPoint(reco::PFTrajectoryPoint::HCALEntrance);
0799 if (atHCAL.isValid()) {
0800 eta_hcal = atHCAL.positionREP().eta();
0801 phi_hcal = atHCAL.positionREP().phi();
0802 }
0803 if (atECAL.isValid()) {
0804 eta_ecal = atECAL.positionREP().eta();
0805 phi_ecal = atECAL.positionREP().phi();
0806 }
0807 }
0808 const auto& ref = ((const reco::PFBlockElementTrack*)&orig)->trackRef();
0809 pt = ref->pt();
0810 px = ref->px();
0811 py = ref->py();
0812 pz = ref->pz();
0813 eta = ref->eta();
0814 phi = ref->phi();
0815 energy = ref->p();
0816 charge = ref->charge();
0817 num_hits = ref->recHitsSize();
0818
0819 reco::MuonRef muonRef = orig.muonRef();
0820 if (muonRef.isNonnull()) {
0821 reco::TrackRef standAloneMu = muonRef->standAloneMuon();
0822 if (standAloneMu.isNonnull()) {
0823 muon_dt_hits = standAloneMu->hitPattern().numberOfValidMuonDTHits();
0824 muon_csc_hits = standAloneMu->hitPattern().numberOfValidMuonCSCHits();
0825 }
0826 muon_type = muonRef->type();
0827 }
0828
0829 } else if (type == reco::PFBlockElement::BREM) {
0830 const auto* orig2 = (const reco::PFBlockElementBrem*)&orig;
0831 const auto& ref = orig2->GsftrackRef();
0832 trajpoint = orig2->indTrajPoint();
0833 if (ref.isNonnull()) {
0834 deltap = orig2->DeltaP();
0835 sigmadeltap = orig2->SigmaDeltaP();
0836 pt = ref->pt();
0837 px = ref->px();
0838 py = ref->py();
0839 pz = ref->pz();
0840 eta = ref->eta();
0841 phi = ref->phi();
0842 energy = ref->p();
0843 charge = ref->charge();
0844 }
0845
0846 const auto& gsfextraref = ref->extra();
0847 if (gsfextraref.isAvailable() && gsfextraref->seedRef().isAvailable()) {
0848 reco::ElectronSeedRef seedref = gsfextraref->seedRef().castTo<reco::ElectronSeedRef>();
0849 if (seedref.isAvailable()) {
0850 if (seedref->isEcalDriven()) {
0851 gsf_electronseed_trkorecal = 1.0;
0852 } else if (seedref->isTrackerDriven()) {
0853 gsf_electronseed_trkorecal = 2.0;
0854 }
0855 }
0856 }
0857
0858 } else if (type == reco::PFBlockElement::GSF) {
0859
0860 const auto* orig2 = (const reco::PFBlockElementGsfTrack*)&orig;
0861 const auto& vec = orig2->Pin();
0862 pt = vec.pt();
0863 px = vec.px();
0864 py = vec.py();
0865 pz = vec.pz();
0866 eta = vec.eta();
0867 phi = vec.phi();
0868 energy = vec.energy();
0869
0870 const auto& vec2 = orig2->Pout();
0871 eta_ecal = vec2.eta();
0872 phi_ecal = vec2.phi();
0873
0874 if (!orig2->GsftrackRefPF().isNull()) {
0875 charge = orig2->GsftrackRefPF()->charge();
0876 num_hits = orig2->GsftrackRefPF()->PFRecBrem().size();
0877 }
0878
0879 const auto& ref = orig2->GsftrackRef();
0880
0881 const auto& gsfextraref = ref->extra();
0882 if (gsfextraref.isAvailable() && gsfextraref->seedRef().isAvailable()) {
0883 reco::ElectronSeedRef seedref = gsfextraref->seedRef().castTo<reco::ElectronSeedRef>();
0884 if (seedref.isAvailable()) {
0885 if (seedref->isEcalDriven()) {
0886 gsf_electronseed_trkorecal = 1.0;
0887 } else if (seedref->isTrackerDriven()) {
0888 gsf_electronseed_trkorecal = 2.0;
0889 }
0890 }
0891 };
0892
0893 } else if (type == reco::PFBlockElement::ECAL || type == reco::PFBlockElement::PS1 ||
0894 type == reco::PFBlockElement::PS2 || type == reco::PFBlockElement::HCAL ||
0895 type == reco::PFBlockElement::HO || type == reco::PFBlockElement::HFHAD ||
0896 type == reco::PFBlockElement::HFEM) {
0897 const auto& ref = ((const reco::PFBlockElementCluster*)&orig)->clusterRef();
0898 if (ref.isNonnull()) {
0899 cluster_flags = ref->flags();
0900 eta = ref->eta();
0901 phi = ref->phi();
0902 pt = ref->pt();
0903 px = ref->position().x();
0904 py = ref->position().y();
0905 pz = ref->position().z();
0906 energy = ref->energy();
0907 corr_energy = ref->correctedEnergy();
0908 layer = ref->layer();
0909 depth = ref->depth();
0910 num_hits = ref->recHitFractions().size();
0911 }
0912 } else if (type == reco::PFBlockElement::SC) {
0913 const auto& clref = ((const reco::PFBlockElementSuperCluster*)&orig)->superClusterRef();
0914 if (clref.isNonnull()) {
0915 cluster_flags = clref->flags();
0916 eta = clref->eta();
0917 phi = clref->phi();
0918 px = clref->position().x();
0919 py = clref->position().y();
0920 pz = clref->position().z();
0921 energy = clref->energy();
0922 num_hits = clref->clustersSize();
0923 }
0924 }
0925 vector<int> tps;
0926 for (const auto& t : trackingparticle_to_element) {
0927 if (t.second == (int)ielem) {
0928 tps.push_back(t.first);
0929 }
0930 }
0931 vector<int> scs;
0932 for (const auto& t : simcluster_to_element) {
0933 if (t.second == (int)ielem) {
0934 scs.push_back(t.first);
0935 }
0936 }
0937
0938 element_pt_.push_back(pt);
0939 element_px_.push_back(px);
0940 element_py_.push_back(py);
0941 element_pz_.push_back(pz);
0942 element_deltap_.push_back(deltap);
0943 element_sigmadeltap_.push_back(sigmadeltap);
0944 element_eta_.push_back(eta);
0945 element_phi_.push_back(phi);
0946 element_energy_.push_back(energy);
0947 element_corr_energy_.push_back(corr_energy);
0948 element_eta_ecal_.push_back(eta_ecal);
0949 element_phi_ecal_.push_back(phi_ecal);
0950 element_eta_hcal_.push_back(eta_hcal);
0951 element_phi_hcal_.push_back(phi_hcal);
0952 element_charge_.push_back(charge);
0953 element_type_.push_back(type);
0954 element_layer_.push_back(layer);
0955 element_depth_.push_back(depth);
0956 element_trajpoint_.push_back(trajpoint);
0957 element_muon_dt_hits_.push_back(muon_dt_hits);
0958 element_muon_csc_hits_.push_back(muon_csc_hits);
0959 element_muon_type_.push_back(muon_type);
0960 element_cluster_flags_.push_back(cluster_flags);
0961 element_gsf_electronseed_trkorecal_.push_back(gsf_electronseed_trkorecal);
0962 element_num_hits_.push_back(num_hits);
0963 }
0964
0965
0966 int icandidate = 0;
0967 for (const auto& cand : pfCandidates) {
0968 pfcandidate_eta_.push_back(cand.eta());
0969 pfcandidate_phi_.push_back(cand.phi());
0970 pfcandidate_pt_.push_back(cand.pt());
0971 pfcandidate_px_.push_back(cand.px());
0972 pfcandidate_py_.push_back(cand.py());
0973 pfcandidate_pz_.push_back(cand.pz());
0974 pfcandidate_energy_.push_back(cand.energy());
0975 pfcandidate_pdgid_.push_back(cand.pdgId());
0976
0977 for (const auto& el : cand.elementsInBlocks()) {
0978 const auto idx_block = el.first.index();
0979 unsigned idx_element_in_block = el.second;
0980
0981 int ielem = -1;
0982 for (const auto& elem_with_index : all_elements) {
0983 ielem += 1;
0984 if (elem_with_index.idx_block == idx_block && elem_with_index.idx_elem == idx_element_in_block) {
0985 break;
0986 }
0987 }
0988 assert(ielem != -1);
0989 element_to_candidate.push_back(make_pair(ielem, icandidate));
0990 }
0991
0992 icandidate += 1;
0993 }
0994
0995 ev_event_ = iEvent.id().event();
0996 ev_lumi_ = iEvent.id().luminosityBlock();
0997 ev_run_ = iEvent.id().run();
0998
0999 t_->Fill();
1000 }
1001
1002 void PFAnalysis::processTrackingParticles(const edm::View<TrackingParticle>& trackingParticles,
1003 edm::Handle<edm::View<TrackingParticle>>& trackingParticlesHandle) {
1004 for (unsigned long ntrackingparticle = 0; ntrackingparticle < trackingParticles.size(); ntrackingparticle++) {
1005 const auto& tp = trackingParticles.at(ntrackingparticle);
1006 edm::RefToBase<TrackingParticle> tpref(trackingParticlesHandle, ntrackingparticle);
1007
1008 math::XYZTLorentzVectorD vtx(0, 0, 0, 0);
1009
1010 if (!tp.decayVertices().empty()) {
1011 vtx = tp.decayVertices().at(0)->position();
1012 }
1013 auto orig_vtx = tp.vertex();
1014
1015
1016 trackingparticle_eta_.push_back(tp.p4().eta());
1017 trackingparticle_phi_.push_back(tp.p4().phi());
1018 trackingparticle_pt_.push_back(tp.p4().pt());
1019 trackingparticle_px_.push_back(tp.p4().px());
1020 trackingparticle_py_.push_back(tp.p4().py());
1021 trackingparticle_pz_.push_back(tp.p4().pz());
1022 trackingparticle_energy_.push_back(tp.p4().energy());
1023 trackingparticle_dvx_.push_back(vtx.x());
1024 trackingparticle_dvy_.push_back(vtx.y());
1025 trackingparticle_dvz_.push_back(vtx.z());
1026 trackingparticle_bx_.push_back(tp.eventId().bunchCrossing());
1027 trackingparticle_ev_.push_back(tp.eventId().event());
1028
1029 trackingparticle_ovx_.push_back(orig_vtx.x());
1030 trackingparticle_ovy_.push_back(orig_vtx.y());
1031 trackingparticle_ovz_.push_back(orig_vtx.z());
1032
1033 trackingparticle_pid_.push_back(tp.pdgId());
1034 }
1035 }
1036
1037
1038 int get_index_triu_vector(int i, int j, int n) {
1039 int k = (n * (n - 1) / 2) - (n - i) * ((n - i) - 1) / 2 + j - i - 1;
1040 return k;
1041 }
1042
1043 pair<int, int> get_triu_vector_index(int k, int n) {
1044 int i = n - 2 - floor(sqrt(-8 * k + 4 * n * (n - 1) - 7) / 2.0 - 0.5);
1045 int j = k + i + 1 - n * (n - 1) / 2 + (n - i) * ((n - i) - 1) / 2;
1046 return make_pair(i, j);
1047 }
1048
1049 pair<vector<ElementWithIndex>, vector<tuple<int, int, float>>> PFAnalysis::processBlocks(
1050 const std::vector<reco::PFBlock>& pfBlocks) {
1051 vector<ElementWithIndex> ret;
1052 vector<tuple<int, int, float>> distances;
1053
1054
1055 int iblock = 0;
1056 for (const auto& block : pfBlocks) {
1057 int ielem = 0;
1058 const auto& linkdata = block.linkData();
1059
1060
1061 for (const auto& link : linkdata) {
1062 const auto vecidx = link.first;
1063 const auto dist = link.second.distance;
1064 const auto& ij = get_triu_vector_index(vecidx, block.elements().size());
1065 auto globalindex_i = ij.first + ret.size();
1066 auto globalindex_j = ij.second + ret.size();
1067 distances.push_back(make_tuple(globalindex_i, globalindex_j, dist));
1068 }
1069
1070 for (const auto& elem : block.elements()) {
1071 ElementWithIndex elem_index(elem, iblock, ielem);
1072 ret.push_back(elem_index);
1073 ielem += 1;
1074 }
1075 iblock += 1;
1076 }
1077 return make_pair(ret, distances);
1078
1079 }
1080
1081 void PFAnalysis::associateClusterToSimCluster(const vector<ElementWithIndex>& all_elements) {
1082 vector<map<uint64_t, double>> detids_elements;
1083 map<uint64_t, double> rechits_energy_all;
1084
1085 int idx_element = 0;
1086 for (const auto& elem : all_elements) {
1087 map<uint64_t, double> detids;
1088 const auto& type = elem.orig.type();
1089
1090 if (type == reco::PFBlockElement::ECAL || type == reco::PFBlockElement::HCAL || type == reco::PFBlockElement::PS1 ||
1091 type == reco::PFBlockElement::PS2 || type == reco::PFBlockElement::HO || type == reco::PFBlockElement::HFHAD ||
1092 type == reco::PFBlockElement::HFEM) {
1093 const auto& clref = elem.orig.clusterRef();
1094 assert(clref.isNonnull());
1095 const auto& cluster = *clref;
1096
1097
1098 const vector<reco::PFRecHitFraction>& rechit_fracs = cluster.recHitFractions();
1099 for (const auto& rh : rechit_fracs) {
1100 const reco::PFRecHit pfrh = *rh.recHitRef();
1101 if (detids.find(pfrh.detId()) != detids.end()) {
1102 continue;
1103 }
1104 detids[pfrh.detId()] += pfrh.energy() * rh.fraction();
1105 const auto id = DetId(pfrh.detId());
1106 float x = 0;
1107 float y = 0;
1108 float z = 0;
1109 float eta = 0;
1110 float phi = 0;
1111
1112 const auto& pos = getHitPosition(id);
1113 x = pos.x();
1114 y = pos.y();
1115 z = pos.z();
1116 eta = pos.eta();
1117 phi = pos.phi();
1118
1119 rechit_x_.push_back(x);
1120 rechit_y_.push_back(y);
1121 rechit_z_.push_back(z);
1122 rechit_det_.push_back(id.det());
1123 rechit_subdet_.push_back(id.subdetId());
1124 rechit_eta_.push_back(eta);
1125 rechit_phi_.push_back(phi);
1126 rechit_e_.push_back(pfrh.energy() * rh.fraction());
1127 rechit_idx_element_.push_back(idx_element);
1128 rechit_detid_.push_back(id.rawId());
1129 rechits_energy_all[id.rawId()] += pfrh.energy() * rh.fraction();
1130 }
1131 } else if (type == reco::PFBlockElement::SC) {
1132 const auto& clref = ((const reco::PFBlockElementSuperCluster*)&(elem.orig))->superClusterRef();
1133 assert(clref.isNonnull());
1134 const auto& cluster = *clref;
1135
1136
1137 const auto& rechit_fracs = cluster.hitsAndFractions();
1138 for (const auto& rh : rechit_fracs) {
1139 if (detids.find(rh.first.rawId()) != detids.end()) {
1140 continue;
1141 }
1142 detids[rh.first.rawId()] += cluster.energy() * rh.second;
1143 const auto id = rh.first;
1144 float x = 0;
1145 float y = 0;
1146 float z = 0;
1147 float eta = 0;
1148 float phi = 0;
1149
1150 const auto& pos = getHitPosition(id);
1151 x = pos.x();
1152 y = pos.y();
1153 z = pos.z();
1154 eta = pos.eta();
1155 phi = pos.phi();
1156
1157 rechit_x_.push_back(x);
1158 rechit_y_.push_back(y);
1159 rechit_z_.push_back(z);
1160 rechit_det_.push_back(id.det());
1161 rechit_subdet_.push_back(id.subdetId());
1162 rechit_eta_.push_back(eta);
1163 rechit_phi_.push_back(phi);
1164 rechit_e_.push_back(rh.second);
1165 rechit_idx_element_.push_back(idx_element);
1166 rechit_detid_.push_back(id.rawId());
1167 rechits_energy_all[id.rawId()] += cluster.energy() * rh.second;
1168 }
1169 }
1170 detids_elements.push_back(detids);
1171 idx_element += 1;
1172 }
1173
1174
1175 int ielement = 0;
1176 for (const auto& detids : detids_elements) {
1177 int isimcluster = 0;
1178 if (!detids.empty()) {
1179 for (const auto& simcluster_detids : simcluster_detids_) {
1180
1181 double cmp = detid_compare(detids, simcluster_detids);
1182 if (cmp > 0) {
1183 simcluster_to_element.push_back(make_pair(isimcluster, ielement));
1184 simcluster_to_element_cmp.push_back((float)cmp);
1185 }
1186 isimcluster += 1;
1187 }
1188 }
1189 ielement += 1;
1190 }
1191 }
1192
1193 void PFAnalysis::beginRun(edm::Run const& iEvent, edm::EventSetup const& es) {
1194 hcons = &es.getData(hcalDDDrecToken_);
1195 aField_ = &es.getData(magFieldToken_);
1196 }
1197
1198 void PFAnalysis::endRun(edm::Run const& iEvent, edm::EventSetup const&) {}
1199
1200 void PFAnalysis::beginJob() { ; }
1201
1202 void PFAnalysis::endJob() {}
1203
1204 void PFAnalysis::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
1205 edm::ParameterSetDescription desc;
1206 desc.setUnknown();
1207 descriptions.addDefault(desc);
1208 }
1209
1210 DEFINE_FWK_MODULE(PFAnalysis);
