[HGCal] Updating HGCalValidator with SimHits, Digis, RecHits and MultiClusters

Validation/HGCalValidation/interface/HGCalValidator.h

@@ -16,6 +16,7 @@
 #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
+#include "DataFormats/ParticleFlowReco/interface/HGCalMultiCluster.h"
 #include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
 #include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
 #include "SimDataFormats/Vertex/interface/SimVertex.h"
@@ -58,12 +59,15 @@ class HGCalValidator : public DQMGlobalEDAnalyzer<HGCalValidatorHistograms> {
   const bool SaveGeneralInfo_;
   const bool doCaloParticlePlots_;
   const bool dolayerclustersPlots_;
+  const bool domulticlustersPlots_;
   const edm::FileInPath cummatbudinxo_;
 
-  std::vector<edm::EDGetTokenT<reco::CaloClusterCollection> > labelToken;
-  edm::EDGetTokenT<std::vector<CaloParticle> > label_cp_effic;
-  edm::EDGetTokenT<std::vector<CaloParticle> > label_cp_fake;
-  edm::EDGetTokenT<std::vector<SimVertex> > simVertices_;
+  std::vector<edm::EDGetTokenT<reco::CaloClusterCollection>> labelToken;
+  edm::EDGetTokenT<reco::CaloClusterCollection> layerclusters_;
+  edm::EDGetTokenT<std::vector<reco::HGCalMultiCluster>> multiClusters_;
+  edm::EDGetTokenT<std::vector<CaloParticle>> label_cp_effic;
+  edm::EDGetTokenT<std::vector<CaloParticle>> label_cp_fake;
+  edm::EDGetTokenT<std::vector<SimVertex>> simVertices_;
   edm::EDGetTokenT<HGCRecHitCollection> recHitsEE_;
   edm::EDGetTokenT<HGCRecHitCollection> recHitsFH_;
   edm::EDGetTokenT<HGCRecHitCollection> recHitsBH_;

Validation/HGCalValidation/interface/HGVHistoProducerAlgo.h

@@ -17,6 +17,7 @@
 #include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
 #include "DataFormats/HGCRecHit/interface/HGCRecHit.h"
+#include "DataFormats/ParticleFlowReco/interface/HGCalMultiCluster.h"
 
 #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
 #include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
@@ -90,6 +91,66 @@ struct HGVHistoProducerAlgoHistograms {
   std::unordered_map<int, ConcurrentMonitorElement> h_caloparticle_energy;
   std::unordered_map<int, ConcurrentMonitorElement> h_caloparticle_pt;
   std::unordered_map<int, ConcurrentMonitorElement> h_caloparticle_phi;
+
+  //For multiclusters
+  ConcurrentMonitorElement h_score_multicl2caloparticle;
+  ConcurrentMonitorElement h_score_caloparticle2multicl;
+  ConcurrentMonitorElement h_energy_vs_score_multicl2caloparticle;
+  ConcurrentMonitorElement h_energy_vs_score_caloparticle2multicl;
+  //contiguous multiclusters
+  ConcurrentMonitorElement h_score_contimulticl2caloparticle;
+  ConcurrentMonitorElement h_score_caloparticle2contimulticl;
+  ConcurrentMonitorElement h_energy_vs_score_contimulticl2caloparticle;
+  ConcurrentMonitorElement h_energy_vs_score_caloparticle2contimulticl;
+  //non contiguous multiclusters
+  ConcurrentMonitorElement h_score_noncontimulticl2caloparticle;
+  ConcurrentMonitorElement h_score_caloparticle2noncontimulticl;
+  ConcurrentMonitorElement h_energy_vs_score_noncontimulticl2caloparticle;
+  ConcurrentMonitorElement h_energy_vs_score_caloparticle2noncontimulticl;
+  //Back to all multiclusters
+  ConcurrentMonitorElement h_num_multicl_eta;
+  ConcurrentMonitorElement h_num_multicl_phi;
+  ConcurrentMonitorElement h_numMerge_multicl_eta;
+  ConcurrentMonitorElement h_numMerge_multicl_phi;
+  ConcurrentMonitorElement h_sharedenergy_multicl2caloparticle;
+  ConcurrentMonitorElement h_sharedenergy_caloparticle2multicl;
+  ConcurrentMonitorElement h_sharedenergy_multicl2caloparticle_vs_eta;
+  ConcurrentMonitorElement h_sharedenergy_multicl2caloparticle_vs_phi;
+  ConcurrentMonitorElement h_sharedenergy_caloparticle2multicl_vs_eta;
+  ConcurrentMonitorElement h_sharedenergy_caloparticle2multicl_vs_phi;
+  ConcurrentMonitorElement h_denom_multicl_eta;
+  ConcurrentMonitorElement h_denom_multicl_phi;
+  ConcurrentMonitorElement h_num_caloparticle_eta;
+  ConcurrentMonitorElement h_num_caloparticle_phi;
+  ConcurrentMonitorElement h_numDup_caloparticle_eta;
+  ConcurrentMonitorElement h_numDup_caloparticle_phi;
+  ConcurrentMonitorElement h_denom_caloparticle_eta;
+  ConcurrentMonitorElement h_denom_caloparticle_phi;
+  ConcurrentMonitorElement h_cellAssociation;
+  ConcurrentMonitorElement h_multiclusternum;
+  ConcurrentMonitorElement h_contmulticlusternum;
+  ConcurrentMonitorElement h_noncontmulticlusternum;
+  ConcurrentMonitorElement h_clusternum_in_multicluster;
+  std::unordered_map<int, ConcurrentMonitorElement> h_clusternum_in_multicluster_perlayer;
+  ConcurrentMonitorElement h_multiplicityOfLCinMCL;
+  ConcurrentMonitorElement h_multiplicity_numberOfEventsHistogram;
+  ConcurrentMonitorElement h_multiplicity_zminus_numberOfEventsHistogram;
+  ConcurrentMonitorElement h_multiplicity_zplus_numberOfEventsHistogram;
+  ConcurrentMonitorElement h_multiplicityOfLCinMCL_vs_layercluster;
+  ConcurrentMonitorElement h_multiplicityOfLCinMCL_vs_layercluster_zminus;
+  ConcurrentMonitorElement h_multiplicityOfLCinMCL_vs_layercluster_zplus;
+  ConcurrentMonitorElement h_multiplicityOfLCinMCL_vs_layerclusterenergy;
+  ConcurrentMonitorElement h_clusternum_in_multicluster_vs_layer;
+  ConcurrentMonitorElement h_multicluster_pt;
+  ConcurrentMonitorElement h_multicluster_eta;
+  ConcurrentMonitorElement h_multicluster_phi;
+  ConcurrentMonitorElement h_multicluster_energy;
+  ConcurrentMonitorElement h_multicluster_x;
+  ConcurrentMonitorElement h_multicluster_y;
+  ConcurrentMonitorElement h_multicluster_z;
+  ConcurrentMonitorElement h_multicluster_firstlayer;
+  ConcurrentMonitorElement h_multicluster_lastlayer;
+  ConcurrentMonitorElement h_multicluster_layersnum;
 };
 
 using Density = hgcal_clustering::Density;
@@ -109,19 +170,24 @@ class HGVHistoProducerAlgo {
                          unsigned layers,
                          std::vector<int> thicknesses,
                          std::string pathtomatbudfile);
+  void bookMultiClusterHistos(DQMStore::ConcurrentBooker& ibook, Histograms& histograms, unsigned layers);
   void layerClusters_to_CaloParticles(const Histograms& histograms,
                                       const reco::CaloClusterCollection& clusters,
                                       std::vector<CaloParticle> const& cP,
                                       std::map<DetId, const HGCRecHit*> const&,
                                       unsigned layers) const;
+  void multiClusters_to_CaloParticles(const Histograms& histograms,
+                                      const std::vector<reco::HGCalMultiCluster>& multiClusters,
+                                      std::vector<CaloParticle> const& cP,
+                                      std::map<DetId, const HGCRecHit*> const&,
+                                      unsigned layers,
+                                      std::vector<bool> contimulti) const;
   void fill_info_histos(const Histograms& histograms, unsigned layers) const;
   void fill_caloparticle_histos(const Histograms& histograms,
                                 int pdgid,
                                 const CaloParticle& caloparticle,
                                 std::vector<SimVertex> const& simVertices) const;
-
   void fill_cluster_histos(const Histograms& histograms, int count, const reco::CaloCluster& cluster) const;
-
   void fill_generic_cluster_histos(const Histograms& histograms,
                                    int count,
                                    const reco::CaloClusterCollection& clusters,
@@ -131,7 +197,12 @@ class HGVHistoProducerAlgo {
                                    std::map<double, double> cummatbudg,
                                    unsigned layers,
                                    std::vector<int> thicknesses) const;
-
+  void fill_multi_cluster_histos(const Histograms& histograms,
+                                 int count,
+                                 const std::vector<reco::HGCalMultiCluster>& multiClusters,
+                                 std::vector<CaloParticle> const& cP,
+                                 std::map<DetId, const HGCRecHit*> const&,
+                                 unsigned layers) const;
   double distance2(const double x1, const double y1, const double x2, const double y2) const;
   double distance(const double x1, const double y1, const double x2, const double y2) const;
 
@@ -145,6 +216,13 @@ class HGVHistoProducerAlgo {
     float fraction;
   };
 
+  struct detIdInfoInMultiCluster {
+    bool operator==(const detIdInfoInMultiCluster& o) const { return multiclusterId == o.multiclusterId; };
+    unsigned int multiclusterId;
+    unsigned int clusterId;
+    float fraction;
+  };
+
   struct caloParticleOnLayer {
     unsigned int caloParticleId;
     float energy = 0;
@@ -183,6 +261,8 @@ class HGVHistoProducerAlgo {
   int nintScore_;
   double minSharedEneFrac_, maxSharedEneFrac_;
   int nintSharedEneFrac_;
+  double minMCLSharedEneFrac_, maxMCLSharedEneFrac_;
+  int nintMCLSharedEneFrac_;
   double minTotNClsperthick_, maxTotNClsperthick_;
   int nintTotNClsperthick_;
   double minTotNcellsperthickperlayer_, maxTotNcellsperthickperlayer_;
@@ -201,6 +281,24 @@ class HGVHistoProducerAlgo {
   int nintClEneperthickperlayer_;
   double minCellsEneDensperthick_, maxCellsEneDensperthick_;
   int nintCellsEneDensperthick_;
+  double minTotNMCLs_, maxTotNMCLs_;
+  int nintTotNMCLs_;
+  double minTotNClsinMCLs_, maxTotNClsinMCLs_;
+  int nintTotNClsinMCLs_;
+  double minTotNClsinMCLsperlayer_, maxTotNClsinMCLsperlayer_;
+  int nintTotNClsinMCLsperlayer_;
+  double minMplofLCs_, maxMplofLCs_;
+  int nintMplofLCs_;
+  double minSizeCLsinMCLs_, maxSizeCLsinMCLs_;
+  int nintSizeCLsinMCLs_;
+  double minClEnepermultiplicity_, maxClEnepermultiplicity_;
+  int nintClEnepermultiplicity_;
+  double minX_, maxX_;
+  int nintX_;
+  double minY_, maxY_;
+  int nintY_;
+  double minZ_, maxZ_;
+  int nintZ_;
 };
 
 #endif

Validation/HGCalValidation/plugins/HGCalDigiValidation.cc

@@ -384,33 +384,37 @@ void HGCalDigiValidation::bookHistograms(DQMStore::IBooker& iB, edm::Run const&,
   std::ostringstream histoname;
   for (int il = 0; il < layers_; ++il) {
     int ilayer = firstLayer_ + il;
+    auto istr1 = std::to_string(ilayer);
+    while (istr1.size() < 2) {
+      istr1.insert(0, "0");
+    }
     histoname.str("");
     histoname << "TOA_"
-              << "layer_" << ilayer;
+              << "layer_" << istr1;
     TOA_.push_back(iB.book1D(histoname.str().c_str(), "toa_", 1024, 0, 1024));
 
     histoname.str("");
     histoname << "ADC_"
-              << "layer_" << ilayer;
+              << "layer_" << istr1;
     ADC_.push_back(iB.book1D(histoname.str().c_str(), "ADCDigiOccupancy", 1024, 0, 1024));
 
     histoname.str("");
     histoname << "TOT_"
-              << "layer_" << ilayer;
+              << "layer_" << istr1;
     TOT_.push_back(iB.book1D(histoname.str().c_str(), "TOTDigiOccupancy", 4096, 0, 4096));
 
     histoname.str("");
     histoname << "DigiOccupancy_XY_"
-              << "layer_" << ilayer;
+              << "layer_" << istr1;
     DigiOccupancy_XY_.push_back(iB.book2D(histoname.str().c_str(), "DigiOccupancy", 50, -500, 500, 50, -500, 500));
 
     histoname.str("");
     histoname << "DigiOccupancy_Plus_"
-              << "layer_" << ilayer;
+              << "layer_" << istr1;
     DigiOccupancy_Plus_.push_back(iB.book1D(histoname.str().c_str(), "DigiOccupancy +z", 100, 0, 1000));
     histoname.str("");
     histoname << "DigiOccupancy_Minus_"
-              << "layer_" << ilayer;
+              << "layer_" << istr1;
     DigiOccupancy_Minus_.push_back(iB.book1D(histoname.str().c_str(), "DigiOccupancy -z", 100, 0, 1000));
   }
 

Validation/HGCalValidation/plugins/HGCalRecHitValidation.cc

@@ -294,25 +294,29 @@ void HGCalRecHitValidation::bookHistograms(DQMStore::IBooker& iB, edm::Run const
   std::ostringstream histoname;
   for (unsigned int il = 0; il < layers_; ++il) {
     int ilayer = firstLayer_ + (int)(il);
+    auto istr1 = std::to_string(ilayer);
+    while (istr1.size() < 2) {
+      istr1.insert(0, "0");
+    }
     histoname.str("");
-    histoname << "HitOccupancy_Plus_layer_" << ilayer;
+    histoname << "HitOccupancy_Plus_layer_" << istr1;
     HitOccupancy_Plus_.push_back(iB.book1D(histoname.str().c_str(), "RecHitOccupancy_Plus", 100, 0, 10000));
     histoname.str("");
-    histoname << "HitOccupancy_Minus_layer_" << ilayer;
+    histoname << "HitOccupancy_Minus_layer_" << istr1;
     HitOccupancy_Minus_.push_back(iB.book1D(histoname.str().c_str(), "RecHitOccupancy_Minus", 100, 0, 10000));
 
     histoname.str("");
     histoname << "EtaPhi_Plus_"
-              << "layer_" << ilayer;
+              << "layer_" << istr1;
     EtaPhi_Plus_.push_back(iB.book2D(histoname.str().c_str(), "Occupancy", 31, 1.45, 3.0, 72, -CLHEP::pi, CLHEP::pi));
     histoname.str("");
     histoname << "EtaPhi_Minus_"
-              << "layer_" << ilayer;
+              << "layer_" << istr1;
     EtaPhi_Minus_.push_back(
         iB.book2D(histoname.str().c_str(), "Occupancy", 31, -3.0, -1.45, 72, -CLHEP::pi, CLHEP::pi));
 
     histoname.str("");
-    histoname << "energy_layer_" << ilayer;
+    histoname << "energy_layer_" << istr1;
     energy_.push_back(iB.book1D(histoname.str().c_str(), "energy_", 500, 0, 1));
   }  //loop over layers ends here
 

Validation/HGCalValidation/plugins/HGCalSimHitValidation.cc

@@ -438,26 +438,30 @@ void HGCalSimHitValidation::bookHistograms(DQMStore::IBooker& iB, edm::Run const
   std::ostringstream histoname;
   for (unsigned int il = 0; il < layers_; ++il) {
     int ilayer = firstLayer_ + static_cast<int>(il);
+    auto istr1 = std::to_string(ilayer);
+    while (istr1.size() < 2) {
+      istr1.insert(0, "0");
+    }
     histoname.str("");
-    histoname << "HitOccupancy_Plus_layer_" << ilayer;
+    histoname << "HitOccupancy_Plus_layer_" << istr1;
     HitOccupancy_Plus_.push_back(iB.book1D(histoname.str().c_str(), "HitOccupancy_Plus", 501, -0.5, 500.5));
     histoname.str("");
-    histoname << "HitOccupancy_Minus_layer_" << ilayer;
+    histoname << "HitOccupancy_Minus_layer_" << istr1;
     HitOccupancy_Minus_.push_back(iB.book1D(histoname.str().c_str(), "HitOccupancy_Minus", 501, -0.5, 500.5));
 
     histoname.str("");
     histoname << "EtaPhi_Plus_"
-              << "layer_" << ilayer;
+              << "layer_" << istr1;
     EtaPhi_Plus_.push_back(iB.book2D(histoname.str().c_str(), "Occupancy", 31, 1.45, 3.0, 72, -CLHEP::pi, CLHEP::pi));
     histoname.str("");
     histoname << "EtaPhi_Minus_"
-              << "layer_" << ilayer;
+              << "layer_" << istr1;
     EtaPhi_Minus_.push_back(
         iB.book2D(histoname.str().c_str(), "Occupancy", 31, -3.0, -1.45, 72, -CLHEP::pi, CLHEP::pi));
 
     for (unsigned int itimeslice = 0; itimeslice < nTimes_; itimeslice++) {
       histoname.str("");
-      histoname << "energy_time_" << itimeslice << "_layer_" << ilayer;
+      histoname << "energy_time_" << itimeslice << "_layer_" << istr1;
       energy_[itimeslice].push_back(iB.book1D(histoname.str().c_str(), "energy_", 100, 0, 0.1));
     }
   }