Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-02-12 09:08:03

0001 // Based on RecoNtuple/HGCalAnalysis with modifications for PF
0002 //
0003 // user include files
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         //rechit energy times simhit fraction
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   // ----------member data ---------------------------
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   // and also the magnetic field
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   // event info
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   //PF Elements
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   //Distance matrix between PF elements
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   //Links between reco, gen and PFCandidate objects
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 }  // constructor
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 }  //clearVariables
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   //Simulated tracks, cleaned up by TrackingTruthAccumulator
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   //Matches reco tracks to sim tracks (TrackingParticle)
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   //Collect all clusters, tracks and superclusters
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   //assert(all_distances.size() > 0);
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   //We need to use the original reco::Track collection for track association
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     //track was not used by PF, we skip as well
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   //Fill genparticles from calorimeter hits
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       //create a map of detId->energy of all the rechits in all the clusters of this SimCluster
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           //compare the two tracks, taking into account that both eventId and trackId need to be compared due to pileup
0746           if (simtrack.eventId() == simtrack2.eventId() && simtrack.trackId() == simtrack2.trackId()) {
0747             simcluster_to_trackingparticle = itp;
0748             //we are satisfied with the first match, in practice there should not be more
0749             break;
0750           }
0751         }  //trackingParticles
0752       }    //simcluster tracks
0753 
0754       simcluster_detids_.push_back(detid_energy);
0755       simcluster_idx_trackingparticle_.push_back(simcluster_to_trackingparticle);
0756 
0757       idx_simcluster += 1;
0758     }  //simclusters
0759   }    //caloParticles
0760 
0761   associateClusterToSimCluster(all_elements);
0762 
0763   //fill elements
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       //requires to keep GsfPFRecTracks
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   //associate candidates to elements
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     }  //elements
0991 
0992     icandidate += 1;
0993   }  //pfCandidates
0994 
0995   ev_event_ = iEvent.id().event();
0996   ev_lumi_ = iEvent.id().luminosityBlock();
0997   ev_run_ = iEvent.id().run();
0998 
0999   t_->Fill();
1000 }  //analyze
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     // fill branches
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 //https://stackoverflow.com/questions/27086195/linear-index-upper-triangular-matrix/27088560
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   //Collect all the elements
1055   int iblock = 0;
1056   for (const auto& block : pfBlocks) {
1057     int ielem = 0;
1058     const auto& linkdata = block.linkData();
1059 
1060     //create a list of global element indices with distances
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     }  //elements
1075     iblock += 1;
1076   }  //blocks
1077   return make_pair(ret, distances);
1078 
1079 }  //processBlocks
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       //all rechits and the energy fractions in this cluster
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       }  //rechit_fracs
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       //all rechits and the energy fractions in this cluster
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       }  //rechit_fracs
1169     }
1170     detids_elements.push_back(detids);
1171     idx_element += 1;
1172   }  //all_elements
1173 
1174   //associate elements (reco clusters) to simclusters
1175   int ielement = 0;
1176   for (const auto& detids : detids_elements) {
1177     int isimcluster = 0;
1178     if (!detids.empty()) {
1179       double sum_e_tot = 0.0;
1180       for (const auto& c : detids) {
1181         sum_e_tot += c.second;
1182       }
1183 
1184       for (const auto& simcluster_detids : simcluster_detids_) {
1185         double sum_e_tot_sc = 0.0;
1186         for (const auto& c : simcluster_detids) {
1187           sum_e_tot_sc += c.second;
1188         }
1189 
1190         //get the energy of the simcluster hits that matches detids of the rechits
1191         double cmp = detid_compare(detids, simcluster_detids);
1192         if (cmp > 0) {
1193           simcluster_to_element.push_back(make_pair(isimcluster, ielement));
1194           simcluster_to_element_cmp.push_back((float)cmp);
1195         }
1196         isimcluster += 1;
1197       }
1198     }  //element had rechits
1199     ielement += 1;
1200   }  //rechit clusters
1201 }
1202 
1203 void PFAnalysis::beginRun(edm::Run const& iEvent, edm::EventSetup const& es) {
1204   hcons = &es.getData(hcalDDDrecToken_);
1205   aField_ = &es.getData(magFieldToken_);
1206 }
1207 
1208 void PFAnalysis::endRun(edm::Run const& iEvent, edm::EventSetup const&) {}
1209 
1210 void PFAnalysis::beginJob() { ; }
1211 
1212 void PFAnalysis::endJob() {}
1213 
1214 void PFAnalysis::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
1215   edm::ParameterSetDescription desc;
1216   desc.setUnknown();
1217   descriptions.addDefault(desc);
1218 }
1219 
1220 DEFINE_FWK_MODULE(PFAnalysis);