Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-07-09 05:00:39 

0001 #include <numeric>
0002 #include <iomanip>
0003 #include <sstream>
0004 
0005 #include "Validation/HGCalValidation/interface/BarrelVHistoProducerAlgo.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
0008 #include "TMath.h"
0009 #include "TLatex.h"
0010 #include "TF1.h"
0011 
0012 using namespace std;
0013 
0014 //Parameters for the score cut. Later, this will become part of the
0015 //configuration parameter for the Barrel associator.
0016 const double ScoreCutLCtoCP_ = 0.1;
0017 const double ScoreCutCPtoLC_ = 0.1;
0018 const double ScoreCutLCtoSC_ = 0.1;
0019 const double ScoreCutSCtoLC_ = 0.1;
0020 
0021 BarrelVHistoProducerAlgo::BarrelVHistoProducerAlgo(const edm::ParameterSet& pset)
0022     :  //parameters for eta
0023       minEta_(pset.getParameter<double>("minEta")),
0024       maxEta_(pset.getParameter<double>("maxEta")),
0025       nintEta_(pset.getParameter<int>("nintEta")),
0026       useFabsEta_(pset.getParameter<bool>("useFabsEta")),
0027 
0028       //parameters for energy
0029       minEne_(pset.getParameter<double>("minEne")),
0030       maxEne_(pset.getParameter<double>("maxEne")),
0031       nintEne_(pset.getParameter<int>("nintEne")),
0032 
0033       //parameters for pt
0034       minPt_(pset.getParameter<double>("minPt")),
0035       maxPt_(pset.getParameter<double>("maxPt")),
0036       nintPt_(pset.getParameter<int>("nintPt")),
0037 
0038       //parameters for phi
0039       minPhi_(pset.getParameter<double>("minPhi")),
0040       maxPhi_(pset.getParameter<double>("maxPhi")),
0041       nintPhi_(pset.getParameter<int>("nintPhi")),
0042 
0043       //parameters for the total amount of energy clustered by all layer clusters (fraction over CaloParticless)
0044       minEneCl_(pset.getParameter<double>("minEneCl")),
0045       maxEneCl_(pset.getParameter<double>("maxEneCl")),
0046       nintEneCl_(pset.getParameter<int>("nintEneCl")),
0047 
0048       //parameters for z positionof vertex plots
0049       minZpos_(pset.getParameter<double>("minZpos")),
0050       maxZpos_(pset.getParameter<double>("maxZpos")),
0051       nintZpos_(pset.getParameter<int>("nintZpos")),
0052 
0053       //Parameters for the total number of SimClusters per layer
0054       minTotNsimClsperlay_(pset.getParameter<double>("minTotNsimClsperlay")),
0055       maxTotNsimClsperlay_(pset.getParameter<double>("maxTotNsimClsperlay")),
0056       nintTotNsimClsperlay_(pset.getParameter<int>("nintTotNsimClsperlay")),
0057 
0058       //Parameters for the total number of layer clusters per layer
0059       minTotNClsperlay_(pset.getParameter<double>("minTotNClsperlay")),
0060       maxTotNClsperlay_(pset.getParameter<double>("maxTotNClsperlay")),
0061       nintTotNClsperlay_(pset.getParameter<int>("nintTotNClsperlay")),
0062 
0063       //Parameters for the energy clustered by layer clusters per layer (fraction over CaloParticless)
0064       minEneClperlay_(pset.getParameter<double>("minEneClperlay")),
0065       maxEneClperlay_(pset.getParameter<double>("maxEneClperlay")),
0066       nintEneClperlay_(pset.getParameter<int>("nintEneClperlay")),
0067 
0068       //Parameters for the score both for:
0069       //1. calo particle to layer clusters association per layer
0070       //2. layer cluster to calo particles association per layer
0071       minScore_(pset.getParameter<double>("minScore")),
0072       maxScore_(pset.getParameter<double>("maxScore")),
0073       nintScore_(pset.getParameter<int>("nintScore")),
0074 
0075       //Parameters for shared energy fraction. That is:
0076       //1. Fraction of each of the layer clusters energy related to a
0077       //calo particle over that calo particle's energy.
0078       //2. Fraction of each of the calo particles energy
0079       //related to a layer cluster over that layer cluster's energy.
0080       minSharedEneFrac_(pset.getParameter<double>("minSharedEneFrac")),
0081       maxSharedEneFrac_(pset.getParameter<double>("maxSharedEneFrac")),
0082       nintSharedEneFrac_(pset.getParameter<int>("nintSharedEneFrac")),
0083 
0084       //parameters for x
0085       minX_(pset.getParameter<double>("minX")),
0086       maxX_(pset.getParameter<double>("maxX")),
0087       nintX_(pset.getParameter<int>("nintX")),
0088 
0089       //parameters for y
0090       minY_(pset.getParameter<double>("minY")),
0091       maxY_(pset.getParameter<double>("maxY")),
0092       nintY_(pset.getParameter<int>("nintY")),
0093 
0094       //parameters for z
0095       minZ_(pset.getParameter<double>("minZ")),
0096       maxZ_(pset.getParameter<double>("maxZ")),
0097       nintZ_(pset.getParameter<int>("nintZ")) {}
0098 
0099 BarrelVHistoProducerAlgo::~BarrelVHistoProducerAlgo() {}
0100 
0101 void BarrelVHistoProducerAlgo::bookInfo(DQMStore::IBooker& ibook, Histograms& histograms) {
0102   histograms.lastLayerEB = ibook.bookInt("lastLayerEB");
0103   histograms.lastLayerHB = ibook.bookInt("lastLayerHB");
0104 }
0105 
0106 void BarrelVHistoProducerAlgo::bookCaloParticleHistos(DQMStore::IBooker& ibook,
0107                                                       Histograms& histograms,
0108                                                       int pdgid,
0109                                                       unsigned int layers) {
0110   histograms.h_caloparticle_eta[pdgid] =
0111       ibook.book1D("N of caloparticle vs eta", "N of caloParticles vs eta", nintEta_, minEta_, maxEta_);
0112   histograms.h_caloparticle_eta_Zorigin[pdgid] =
0113       ibook.book2D("Eta vs Zorigin", "Eta vs Zorigin", nintEta_, minEta_, maxEta_, nintZpos_, minZpos_, maxZpos_);
0114 
0115   histograms.h_caloparticle_energy[pdgid] =
0116       ibook.book1D("Energy", "Energy of CaloParticles; Energy [GeV]", nintEne_, minEne_, maxEne_);
0117   histograms.h_caloparticle_pt[pdgid] = ibook.book1D("Pt", "Pt of CaloParticles", nintPt_, minPt_, maxPt_);
0118   histograms.h_caloparticle_phi[pdgid] = ibook.book1D("Phi", "Phi of CaloParticles", nintPhi_, minPhi_, maxPhi_);
0119   histograms.h_caloparticle_selfenergy[pdgid] =
0120       ibook.book1D("SelfEnergy", "Total Energy of Hits in Sim Clusters (matched)", nintEne_, minEne_, maxEne_);
0121   histograms.h_caloparticle_energyDifference[pdgid] =
0122       ibook.book1D("EnergyDifference", "(Energy-SelfEnergy)/Energy", 300., -5., 1.);
0123 
0124   histograms.h_caloparticle_nSimClusters[pdgid] =
0125       ibook.book1D("Num Sim Clusters", "Num Sim Clusters in CaloParticles", 100, 0., 100.);
0126   histograms.h_caloparticle_nHitsInSimClusters[pdgid] =
0127       ibook.book1D("Num Hits in Sim Clusters", "Num Hits in Sim Clusters in CaloParticles", 1000, 0., 1000.);
0128   histograms.h_caloparticle_nHitsInSimClusters_matchedtoRecHit[pdgid] = ibook.book1D(
0129       "Num Rec-matched Hits in Sim Clusters", "Num Hits in Sim Clusters (matched) in CaloParticles", 1000, 0., 1000.);
0130 
0131   histograms.h_caloparticle_nHits_matched_energy[pdgid] =
0132       ibook.book1D("Energy of Rec-matched Hits", "Energy of Hits in Sim Clusters (matched)", 100, 0., 10.);
0133   histograms.h_caloparticle_nHits_matched_energy_layer[pdgid] =
0134       ibook.book2D("Energy of Rec-matched Hits vs layer",
0135                    "Energy of Hits in Sim Clusters (matched) vs layer",
0136                    layers,
0137                    0.,
0138                    (float)layers,
0139                    100,
0140                    0.,
0141                    10.);
0142   histograms.h_caloparticle_nHits_matched_energy_layer_1SimCl[pdgid] =
0143       ibook.book2D("Energy of Rec-matched Hits vs layer (1SC)",
0144                    "Energy of Hits only 1 Sim Clusters (matched) vs layer",
0145                    layers,
0146                    0.,
0147                    (float)layers,
0148                    100,
0149                    0.,
0150                    10.);
0151   histograms.h_caloparticle_sum_energy_layer[pdgid] =
0152       ibook.book2D("Rec-matched Hits Sum Energy vs layer",
0153                    "Rescaled Sum Energy of Hits in Sim Clusters (matched) vs layer",
0154                    layers,
0155                    0.,
0156                    (float)layers,
0157                    110,
0158                    0.,
0159                    110.);
0160   histograms.h_caloparticle_fractions[pdgid] =
0161       ibook.book2D("HitFractions", "Hit fractions;Hit fraction;E_{hit}^{2} fraction", 101, 0, 1.01, 100, 0, 1);
0162   histograms.h_caloparticle_fractions_weight[pdgid] = ibook.book2D(
0163       "HitFractions_weighted", "Hit fractions weighted;Hit fraction;E_{hit}^{2} fraction", 101, 0, 1.01, 100, 0, 1);
0164 
0165   histograms.h_caloparticle_firstlayer[pdgid] =
0166       ibook.book1D("First Layer", "First layer of the CaloParticles", layers, 0., (float)layers);
0167   histograms.h_caloparticle_lastlayer[pdgid] =
0168       ibook.book1D("Last Layer", "Last layer of the CaloParticles", layers, 0., (float)layers);
0169   histograms.h_caloparticle_layersnum[pdgid] =
0170       ibook.book1D("Number of Layers", "Number of layers of the CaloParticles", layers, 0., (float)layers);
0171   histograms.h_caloparticle_firstlayer_matchedtoRecHit[pdgid] = ibook.book1D(
0172       "First Layer (rec-matched hit)", "First layer of the CaloParticles (matched)", layers, 0., (float)layers);
0173   histograms.h_caloparticle_lastlayer_matchedtoRecHit[pdgid] = ibook.book1D(
0174       "Last Layer (rec-matched hit)", "Last layer of the CaloParticles (matched)", layers, 0., (float)layers);
0175   histograms.h_caloparticle_layersnum_matchedtoRecHit[pdgid] =
0176       ibook.book1D("Number of Layers (rec-matched hit)",
0177                    "Number of layers of the CaloParticles (matched)",
0178                    layers,
0179                    0.,
0180                    (float)layers);
0181 }
0182 
0183 void BarrelVHistoProducerAlgo::bookSimClusterHistos(DQMStore::IBooker& ibook,
0184                                                     Histograms& histograms,
0185                                                     unsigned int layers) {
0186   //---------------------------------------------------------------------------------------------------------------------------
0187   for (unsigned ilayer = 0; ilayer < layers; ++ilayer) {
0188     auto istr1 = std::to_string(ilayer);
0189     std::string istr2 = istr1 + " in barrel";
0190     istr1 = "_barrel" + istr1;
0191     histograms.h_simclusternum_perlayer[ilayer] = ibook.book1D("totsimclusternum_layer_" + istr1,
0192                                                                "total number of SimClusters for layer " + istr1,
0193                                                                nintTotNsimClsperlay_,
0194                                                                minTotNsimClsperlay_,
0195                                                                maxTotNsimClsperlay_);
0196 
0197   }  //end of loop over layers
0198 }
0199 
0200 void BarrelVHistoProducerAlgo::bookSimClusterAssociationHistos(DQMStore::IBooker& ibook,
0201                                                                Histograms& histograms,
0202                                                                unsigned int layers) {
0203   std::unordered_map<int, dqm::reco::MonitorElement*> denom_layercl_in_simcl_eta_perlayer;
0204   denom_layercl_in_simcl_eta_perlayer.clear();
0205   std::unordered_map<int, dqm::reco::MonitorElement*> denom_layercl_in_simcl_phi_perlayer;
0206   denom_layercl_in_simcl_phi_perlayer.clear();
0207   std::unordered_map<int, dqm::reco::MonitorElement*> score_layercl2simcluster_perlayer;
0208   score_layercl2simcluster_perlayer.clear();
0209   std::unordered_map<int, dqm::reco::MonitorElement*> sharedenergy_layercl2simcluster_perlayer;
0210   sharedenergy_layercl2simcluster_perlayer.clear();
0211   std::unordered_map<int, dqm::reco::MonitorElement*> energy_vs_score_layercl2simcluster_perlayer;
0212   energy_vs_score_layercl2simcluster_perlayer.clear();
0213   std::unordered_map<int, dqm::reco::MonitorElement*> num_layercl_in_simcl_eta_perlayer;
0214   num_layercl_in_simcl_eta_perlayer.clear();
0215   std::unordered_map<int, dqm::reco::MonitorElement*> num_layercl_in_simcl_phi_perlayer;
0216   num_layercl_in_simcl_phi_perlayer.clear();
0217   std::unordered_map<int, dqm::reco::MonitorElement*> numMerge_layercl_in_simcl_eta_perlayer;
0218   numMerge_layercl_in_simcl_eta_perlayer.clear();
0219   std::unordered_map<int, dqm::reco::MonitorElement*> numMerge_layercl_in_simcl_phi_perlayer;
0220   numMerge_layercl_in_simcl_phi_perlayer.clear();
0221   std::unordered_map<int, dqm::reco::MonitorElement*> sharedenergy_layercl2simcluster_vs_eta_perlayer;
0222   sharedenergy_layercl2simcluster_vs_eta_perlayer.clear();
0223   std::unordered_map<int, dqm::reco::MonitorElement*> sharedenergy_layercl2simcluster_vs_phi_perlayer;
0224   sharedenergy_layercl2simcluster_vs_phi_perlayer.clear();
0225   std::unordered_map<int, dqm::reco::MonitorElement*> denom_simcluster_eta_perlayer;
0226   denom_simcluster_eta_perlayer.clear();
0227   std::unordered_map<int, dqm::reco::MonitorElement*> denom_simcluster_phi_perlayer;
0228   denom_simcluster_phi_perlayer.clear();
0229   std::unordered_map<int, dqm::reco::MonitorElement*> score_simcluster2layercl_perlayer;
0230   score_simcluster2layercl_perlayer.clear();
0231   std::unordered_map<int, dqm::reco::MonitorElement*> sharedenergy_simcluster2layercl_perlayer;
0232   sharedenergy_simcluster2layercl_perlayer.clear();
0233   std::unordered_map<int, dqm::reco::MonitorElement*> energy_vs_score_simcluster2layercl_perlayer;
0234   energy_vs_score_simcluster2layercl_perlayer.clear();
0235   std::unordered_map<int, dqm::reco::MonitorElement*> num_simcluster_eta_perlayer;
0236   num_simcluster_eta_perlayer.clear();
0237   std::unordered_map<int, dqm::reco::MonitorElement*> num_simcluster_phi_perlayer;
0238   num_simcluster_phi_perlayer.clear();
0239   std::unordered_map<int, dqm::reco::MonitorElement*> numDup_simcluster_eta_perlayer;
0240   numDup_simcluster_eta_perlayer.clear();
0241   std::unordered_map<int, dqm::reco::MonitorElement*> numDup_simcluster_phi_perlayer;
0242   numDup_simcluster_phi_perlayer.clear();
0243   std::unordered_map<int, dqm::reco::MonitorElement*> sharedenergy_simcluster2layercl_vs_eta_perlayer;
0244   sharedenergy_simcluster2layercl_vs_eta_perlayer.clear();
0245   std::unordered_map<int, dqm::reco::MonitorElement*> sharedenergy_simcluster2layercl_vs_phi_perlayer;
0246   sharedenergy_simcluster2layercl_vs_phi_perlayer.clear();
0247 
0248   //---------------------------------------------------------------------------------------------------------------------------
0249   for (unsigned ilayer = 0; ilayer < layers; ++ilayer) {
0250     auto istr1 = std::to_string(ilayer);
0251     std::string istr2 = istr1 + " in barrel";
0252     istr1 = "_barrel" + istr1;
0253     //-------------------------------------------------------------------------------------------------------------------------
0254     denom_layercl_in_simcl_eta_perlayer[ilayer] =
0255         ibook.book1D("Denom_LayerCluster_in_SimCluster_Eta_perlayer" + istr1,
0256                      "Denom LayerCluster in SimCluster Eta per Layer Cluster for layer " + istr2,
0257                      nintEta_,
0258                      minEta_,
0259                      maxEta_);
0260     //-------------------------------------------------------------------------------------------------------------------------
0261     denom_layercl_in_simcl_phi_perlayer[ilayer] =
0262         ibook.book1D("Denom_LayerCluster_in_SimCluster_Phi_perlayer" + istr1,
0263                      "Denom LayerCluster in SimCluster Phi per Layer Cluster for layer " + istr2,
0264                      nintPhi_,
0265                      minPhi_,
0266                      maxPhi_);
0267     //-------------------------------------------------------------------------------------------------------------------------
0268     score_layercl2simcluster_perlayer[ilayer] = ibook.book1D("Score_layercl2simcluster_perlayer" + istr1,
0269                                                              "Score of Layer Cluster per SimCluster for layer " + istr2,
0270                                                              nintScore_,
0271                                                              minScore_,
0272                                                              maxScore_);
0273     //-------------------------------------------------------------------------------------------------------------------------
0274     score_simcluster2layercl_perlayer[ilayer] = ibook.book1D("Score_simcluster2layercl_perlayer" + istr1,
0275                                                              "Score of SimCluster per Layer Cluster for layer " + istr2,
0276                                                              nintScore_,
0277                                                              minScore_,
0278                                                              maxScore_);
0279     //-------------------------------------------------------------------------------------------------------------------------
0280     energy_vs_score_simcluster2layercl_perlayer[ilayer] =
0281         ibook.book2D("Energy_vs_Score_simcluster2layer_perlayer" + istr1,
0282                      "Energy vs Score of SimCluster per Layer Cluster for layer " + istr2,
0283                      nintScore_,
0284                      minScore_,
0285                      maxScore_,
0286                      nintSharedEneFrac_,
0287                      minSharedEneFrac_,
0288                      maxSharedEneFrac_);
0289     //-------------------------------------------------------------------------------------------------------------------------
0290     energy_vs_score_layercl2simcluster_perlayer[ilayer] =
0291         ibook.book2D("Energy_vs_Score_layer2simcluster_perlayer" + istr1,
0292                      "Energy vs Score of Layer Cluster per SimCluster for layer " + istr2,
0293                      nintScore_,
0294                      minScore_,
0295                      maxScore_,
0296                      nintSharedEneFrac_,
0297                      minSharedEneFrac_,
0298                      maxSharedEneFrac_);
0299     //-------------------------------------------------------------------------------------------------------------------------
0300     sharedenergy_simcluster2layercl_perlayer[ilayer] =
0301         ibook.book1D("SharedEnergy_simcluster2layercl_perlayer" + istr1,
0302                      "Shared Energy of SimCluster per Layer Cluster for layer " + istr2,
0303                      nintSharedEneFrac_,
0304                      minSharedEneFrac_,
0305                      maxSharedEneFrac_);
0306     //-------------------------------------------------------------------------------------------------------------------------
0307     sharedenergy_simcluster2layercl_vs_eta_perlayer[ilayer] =
0308         ibook.bookProfile("SharedEnergy_simcluster2layercl_vs_eta_perlayer" + istr1,
0309                           "Shared Energy of SimCluster vs #eta per best Layer Cluster for layer " + istr2,
0310                           nintEta_,
0311                           minEta_,
0312                           maxEta_,
0313                           minSharedEneFrac_,
0314                           maxSharedEneFrac_);
0315     //-------------------------------------------------------------------------------------------------------------------------
0316     sharedenergy_simcluster2layercl_vs_phi_perlayer[ilayer] =
0317         ibook.bookProfile("SharedEnergy_simcluster2layercl_vs_phi_perlayer" + istr1,
0318                           "Shared Energy of SimCluster vs #phi per best Layer Cluster for layer " + istr2,
0319                           nintPhi_,
0320                           minPhi_,
0321                           maxPhi_,
0322                           minSharedEneFrac_,
0323                           maxSharedEneFrac_);
0324     //-------------------------------------------------------------------------------------------------------------------------
0325     sharedenergy_layercl2simcluster_perlayer[ilayer] =
0326         ibook.book1D("SharedEnergy_layercluster2simcluster_perlayer" + istr1,
0327                      "Shared Energy of Layer Cluster per SimCluster for layer " + istr2,
0328                      nintSharedEneFrac_,
0329                      minSharedEneFrac_,
0330                      maxSharedEneFrac_);
0331     //-------------------------------------------------------------------------------------------------------------------------
0332     sharedenergy_layercl2simcluster_vs_eta_perlayer[ilayer] =
0333         ibook.bookProfile("SharedEnergy_layercl2simcluster_vs_eta_perlayer" + istr1,
0334                           "Shared Energy of LayerCluster vs #eta per best SimCluster for layer " + istr2,
0335                           nintEta_,
0336                           minEta_,
0337                           maxEta_,
0338                           minSharedEneFrac_,
0339                           maxSharedEneFrac_);
0340     //-------------------------------------------------------------------------------------------------------------------------
0341     sharedenergy_layercl2simcluster_vs_phi_perlayer[ilayer] =
0342         ibook.bookProfile("SharedEnergy_layercl2simcluster_vs_phi_perlayer" + istr1,
0343                           "Shared Energy of LayerCluster vs #phi per best SimCluster for layer " + istr2,
0344                           nintPhi_,
0345                           minPhi_,
0346                           maxPhi_,
0347                           minSharedEneFrac_,
0348                           maxSharedEneFrac_);
0349     //-------------------------------------------------------------------------------------------------------------------------
0350     num_simcluster_eta_perlayer[ilayer] = ibook.book1D("Num_SimCluster_Eta_perlayer" + istr1,
0351                                                        "Num SimCluster Eta per Layer Cluster for layer " + istr2,
0352                                                        nintEta_,
0353                                                        minEta_,
0354                                                        maxEta_);
0355     //-------------------------------------------------------------------------------------------------------------------------
0356     numDup_simcluster_eta_perlayer[ilayer] =
0357         ibook.book1D("NumDup_SimCluster_Eta_perlayer" + istr1,
0358                      "Num Duplicate SimCluster Eta per Layer Cluster for layer " + istr2,
0359                      nintEta_,
0360                      minEta_,
0361                      maxEta_);
0362     //-------------------------------------------------------------------------------------------------------------------------
0363     denom_simcluster_eta_perlayer[ilayer] = ibook.book1D("Denom_SimCluster_Eta_perlayer" + istr1,
0364                                                          "Denom SimCluster Eta per Layer Cluster for layer " + istr2,
0365                                                          nintEta_,
0366                                                          minEta_,
0367                                                          maxEta_);
0368     //-------------------------------------------------------------------------------------------------------------------------
0369     num_simcluster_phi_perlayer[ilayer] = ibook.book1D("Num_SimCluster_Phi_perlayer" + istr1,
0370                                                        "Num SimCluster Phi per Layer Cluster for layer " + istr2,
0371                                                        nintPhi_,
0372                                                        minPhi_,
0373                                                        maxPhi_);
0374     //-------------------------------------------------------------------------------------------------------------------------
0375     numDup_simcluster_phi_perlayer[ilayer] =
0376         ibook.book1D("NumDup_SimCluster_Phi_perlayer" + istr1,
0377                      "Num Duplicate SimCluster Phi per Layer Cluster for layer " + istr2,
0378                      nintPhi_,
0379                      minPhi_,
0380                      maxPhi_);
0381     //-------------------------------------------------------------------------------------------------------------------------
0382     denom_simcluster_phi_perlayer[ilayer] = ibook.book1D("Denom_SimCluster_Phi_perlayer" + istr1,
0383                                                          "Denom SimCluster Phi per Layer Cluster for layer " + istr2,
0384                                                          nintPhi_,
0385                                                          minPhi_,
0386                                                          maxPhi_);
0387     //-------------------------------------------------------------------------------------------------------------------------
0388     num_layercl_in_simcl_eta_perlayer[ilayer] =
0389         ibook.book1D("Num_LayerCluster_in_SimCluster_Eta_perlayer" + istr1,
0390                      "Num LayerCluster Eta per Layer Cluster in SimCluster for layer " + istr2,
0391                      nintEta_,
0392                      minEta_,
0393                      maxEta_);
0394     //-------------------------------------------------------------------------------------------------------------------------
0395     numMerge_layercl_in_simcl_eta_perlayer[ilayer] =
0396         ibook.book1D("NumMerge_LayerCluster_in_SimCluster_Eta_perlayer" + istr1,
0397                      "Num Merge LayerCluster Eta per Layer Cluster in SimCluster for layer " + istr2,
0398                      nintEta_,
0399                      minEta_,
0400                      maxEta_);
0401     //-------------------------------------------------------------------------------------------------------------------------
0402     num_layercl_in_simcl_phi_perlayer[ilayer] =
0403         ibook.book1D("Num_LayerCluster_in_SimCluster_Phi_perlayer" + istr1,
0404                      "Num LayerCluster Phi per Layer Cluster in SimCluster for layer " + istr2,
0405                      nintPhi_,
0406                      minPhi_,
0407                      maxPhi_);
0408     //-------------------------------------------------------------------------------------------------------------------------
0409     numMerge_layercl_in_simcl_phi_perlayer[ilayer] =
0410         ibook.book1D("NumMerge_LayerCluster_in_SimCluster_Phi_perlayer" + istr1,
0411                      "Num Merge LayerCluster Phi per Layer Cluster in SimCluster for layer " + istr2,
0412                      nintPhi_,
0413                      minPhi_,
0414                      maxPhi_);
0415 
0416   }  //end of loop over layers
0417 
0418   histograms.h_denom_layercl_in_simcl_eta_perlayer.push_back(std::move(denom_layercl_in_simcl_eta_perlayer));
0419   histograms.h_denom_layercl_in_simcl_phi_perlayer.push_back(std::move(denom_layercl_in_simcl_phi_perlayer));
0420   histograms.h_score_layercl2simcluster_perlayer.push_back(std::move(score_layercl2simcluster_perlayer));
0421   histograms.h_sharedenergy_layercl2simcluster_perlayer.push_back(std::move(sharedenergy_layercl2simcluster_perlayer));
0422   histograms.h_energy_vs_score_layercl2simcluster_perlayer.push_back(
0423       std::move(energy_vs_score_layercl2simcluster_perlayer));
0424   histograms.h_num_layercl_in_simcl_eta_perlayer.push_back(std::move(num_layercl_in_simcl_eta_perlayer));
0425   histograms.h_num_layercl_in_simcl_phi_perlayer.push_back(std::move(num_layercl_in_simcl_phi_perlayer));
0426   histograms.h_numMerge_layercl_in_simcl_eta_perlayer.push_back(std::move(numMerge_layercl_in_simcl_eta_perlayer));
0427   histograms.h_numMerge_layercl_in_simcl_phi_perlayer.push_back(std::move(numMerge_layercl_in_simcl_phi_perlayer));
0428   histograms.h_sharedenergy_layercl2simcluster_vs_eta_perlayer.push_back(
0429       std::move(sharedenergy_layercl2simcluster_vs_eta_perlayer));
0430   histograms.h_sharedenergy_layercl2simcluster_vs_phi_perlayer.push_back(
0431       std::move(sharedenergy_layercl2simcluster_vs_phi_perlayer));
0432   histograms.h_denom_simcluster_eta_perlayer.push_back(std::move(denom_simcluster_eta_perlayer));
0433   histograms.h_denom_simcluster_phi_perlayer.push_back(std::move(denom_simcluster_phi_perlayer));
0434   histograms.h_score_simcluster2layercl_perlayer.push_back(std::move(score_simcluster2layercl_perlayer));
0435   histograms.h_sharedenergy_simcluster2layercl_perlayer.push_back(std::move(sharedenergy_simcluster2layercl_perlayer));
0436   histograms.h_energy_vs_score_simcluster2layercl_perlayer.push_back(
0437       std::move(energy_vs_score_simcluster2layercl_perlayer));
0438   histograms.h_num_simcluster_eta_perlayer.push_back(std::move(num_simcluster_eta_perlayer));
0439   histograms.h_num_simcluster_phi_perlayer.push_back(std::move(num_simcluster_phi_perlayer));
0440   histograms.h_numDup_simcluster_eta_perlayer.push_back(std::move(numDup_simcluster_eta_perlayer));
0441   histograms.h_numDup_simcluster_phi_perlayer.push_back(std::move(numDup_simcluster_phi_perlayer));
0442   histograms.h_sharedenergy_simcluster2layercl_vs_eta_perlayer.push_back(
0443       std::move(sharedenergy_simcluster2layercl_vs_eta_perlayer));
0444   histograms.h_sharedenergy_simcluster2layercl_vs_phi_perlayer.push_back(
0445       std::move(sharedenergy_simcluster2layercl_vs_phi_perlayer));
0446 }
0447 void BarrelVHistoProducerAlgo::bookClusterHistos_ClusterLevel(DQMStore::IBooker& ibook,
0448                                                               Histograms& histograms,
0449                                                               unsigned int layers) {
0450   //---------------------------------------------------------------------------------------------------------------------------
0451   histograms.h_cluster_eta.push_back(
0452       ibook.book1D("num_reco_cluster_eta", "N of reco clusters vs eta", nintEta_, minEta_, maxEta_));
0453 
0454   //---------------------------------------------------------------------------------------------------------------------------
0455   for (unsigned ilayer = 0; ilayer < layers; ++ilayer) {
0456     auto istr1 = std::to_string(ilayer);
0457     std::string istr2 = istr1 + " in barrel";
0458     istr1 = "_barrel" + istr1;
0459 
0460     histograms.h_clusternum_perlayer[ilayer] = ibook.book1D("totclusternum_layer_" + istr1,
0461                                                             "total number of layer clusters for layer " + istr2,
0462                                                             nintTotNClsperlay_,
0463                                                             minTotNClsperlay_,
0464                                                             maxTotNClsperlay_);
0465     histograms.h_energyclustered_perlayer[ilayer] = ibook.book1D(
0466         "energyclustered_perlayer" + istr1,
0467         "percent of total energy clustered by layer clusters over CaloParticless energy for layer " + istr2,
0468         nintEneClperlay_,
0469         minEneClperlay_,
0470         maxEneClperlay_);
0471   }
0472 }
0473 
0474 void BarrelVHistoProducerAlgo::bookClusterHistos_LCtoCP_association(DQMStore::IBooker& ibook,
0475                                                                     Histograms& histograms,
0476                                                                     unsigned int layers) {
0477   //----------------------------------------------------------------------------------------------------------------------------
0478   for (unsigned ilayer = 0; ilayer < layers; ++ilayer) {
0479     auto istr1 = std::to_string(ilayer);
0480     std::string istr2 = istr1 + " in barrel";
0481     istr1 = "_barrel" + istr1;
0482     histograms.h_score_layercl2caloparticle_perlayer[ilayer] =
0483         ibook.book1D("Score_layercl2caloparticle_perlayer" + istr1,
0484                      "Score of Layer Cluster per CaloParticle for layer " + istr2,
0485                      nintScore_,
0486                      minScore_,
0487                      maxScore_);
0488     histograms.h_score_caloparticle2layercl_perlayer[ilayer] =
0489         ibook.book1D("Score_caloparticle2layercl_perlayer" + istr1,
0490                      "Score of CaloParticle per Layer Cluster for layer " + istr2,
0491                      nintScore_,
0492                      minScore_,
0493                      maxScore_);
0494     histograms.h_energy_vs_score_caloparticle2layercl_perlayer[ilayer] =
0495         ibook.book2D("Energy_vs_Score_caloparticle2layer_perlayer" + istr1,
0496                      "Energy vs Score of CaloParticle per Layer Cluster for layer " + istr2,
0497                      nintScore_,
0498                      minScore_,
0499                      maxScore_,
0500                      nintSharedEneFrac_,
0501                      minSharedEneFrac_,
0502                      maxSharedEneFrac_);
0503     histograms.h_energy_vs_score_layercl2caloparticle_perlayer[ilayer] =
0504         ibook.book2D("Energy_vs_Score_layer2caloparticle_perlayer" + istr1,
0505                      "Energy vs Score of Layer Cluster per CaloParticle Layer for layer " + istr2,
0506                      nintScore_,
0507                      minScore_,
0508                      maxScore_,
0509                      nintSharedEneFrac_,
0510                      minSharedEneFrac_,
0511                      maxSharedEneFrac_);
0512     histograms.h_sharedenergy_caloparticle2layercl_perlayer[ilayer] =
0513         ibook.book1D("SharedEnergy_caloparticle2layercl_perlayer" + istr1,
0514                      "Shared Energy of CaloParticle per Layer Cluster for layer " + istr2,
0515                      nintSharedEneFrac_,
0516                      minSharedEneFrac_,
0517                      maxSharedEneFrac_);
0518     histograms.h_sharedenergy_caloparticle2layercl_vs_eta_perlayer[ilayer] =
0519         ibook.bookProfile("SharedEnergy_caloparticle2layercl_vs_eta_perlayer" + istr1,
0520                           "Shared Energy of CaloParticle vs #eta per best Layer Cluster for layer " + istr2,
0521                           nintEta_,
0522                           minEta_,
0523                           maxEta_,
0524                           minSharedEneFrac_,
0525                           maxSharedEneFrac_);
0526     histograms.h_sharedenergy_caloparticle2layercl_vs_phi_perlayer[ilayer] =
0527         ibook.bookProfile("SharedEnergy_caloparticle2layercl_vs_phi_perlayer" + istr1,
0528                           "Shared Energy of CaloParticle vs #phi per best Layer Cluster for layer " + istr2,
0529                           nintPhi_,
0530                           minPhi_,
0531                           maxPhi_,
0532                           minSharedEneFrac_,
0533                           maxSharedEneFrac_);
0534     histograms.h_sharedenergy_layercl2caloparticle_perlayer[ilayer] =
0535         ibook.book1D("SharedEnergy_layercluster2caloparticle_perlayer" + istr1,
0536                      "Shared Energy of Layer Cluster per Layer Calo Particle for layer " + istr2,
0537                      nintSharedEneFrac_,
0538                      minSharedEneFrac_,
0539                      maxSharedEneFrac_);
0540     histograms.h_sharedenergy_layercl2caloparticle_vs_eta_perlayer[ilayer] =
0541         ibook.bookProfile("SharedEnergy_layercl2caloparticle_vs_eta_perlayer" + istr1,
0542                           "Shared Energy of LayerCluster vs #eta per best Calo Particle for layer " + istr2,
0543                           nintEta_,
0544                           minEta_,
0545                           maxEta_,
0546                           minSharedEneFrac_,
0547                           maxSharedEneFrac_);
0548     histograms.h_sharedenergy_layercl2caloparticle_vs_phi_perlayer[ilayer] =
0549         ibook.bookProfile("SharedEnergy_layercl2caloparticle_vs_phi_perlayer" + istr1,
0550                           "Shared Energy of LayerCluster vs #phi per best Calo Particle for layer " + istr2,
0551                           nintPhi_,
0552                           minPhi_,
0553                           maxPhi_,
0554                           minSharedEneFrac_,
0555                           maxSharedEneFrac_);
0556     histograms.h_num_caloparticle_eta_perlayer[ilayer] =
0557         ibook.book1D("Num_CaloParticle_Eta_perlayer" + istr1,
0558                      "Num CaloParticle Eta per Layer Cluster for layer " + istr2,
0559                      nintEta_,
0560                      minEta_,
0561                      maxEta_);
0562     histograms.h_numDup_caloparticle_eta_perlayer[ilayer] =
0563         ibook.book1D("NumDup_CaloParticle_Eta_perlayer" + istr1,
0564                      "Num Duplicate CaloParticle Eta per Layer Cluster for layer " + istr2,
0565                      nintEta_,
0566                      minEta_,
0567                      maxEta_);
0568     histograms.h_denom_caloparticle_eta_perlayer[ilayer] =
0569         ibook.book1D("Denom_CaloParticle_Eta_perlayer" + istr1,
0570                      "Denom CaloParticle Eta per Layer Cluster for layer " + istr2,
0571                      nintEta_,
0572                      minEta_,
0573                      maxEta_);
0574     histograms.h_num_caloparticle_phi_perlayer[ilayer] =
0575         ibook.book1D("Num_CaloParticle_Phi_perlayer" + istr1,
0576                      "Num CaloParticle Phi per Layer Cluster for layer " + istr2,
0577                      nintPhi_,
0578                      minPhi_,
0579                      maxPhi_);
0580     histograms.h_numDup_caloparticle_phi_perlayer[ilayer] =
0581         ibook.book1D("NumDup_CaloParticle_Phi_perlayer" + istr1,
0582                      "Num Duplicate CaloParticle Phi per Layer Cluster for layer " + istr2,
0583                      nintPhi_,
0584                      minPhi_,
0585                      maxPhi_);
0586     histograms.h_denom_caloparticle_phi_perlayer[ilayer] =
0587         ibook.book1D("Denom_CaloParticle_Phi_perlayer" + istr1,
0588                      "Denom CaloParticle Phi per Layer Cluster for layer " + istr2,
0589                      nintPhi_,
0590                      minPhi_,
0591                      maxPhi_);
0592     histograms.h_num_layercl_eta_perlayer[ilayer] =
0593         ibook.book1D("Num_LayerCluster_Eta_perlayer" + istr1,
0594                      "Num LayerCluster Eta per Layer Cluster for layer " + istr2,
0595                      nintEta_,
0596                      minEta_,
0597                      maxEta_);
0598     histograms.h_numMerge_layercl_eta_perlayer[ilayer] =
0599         ibook.book1D("NumMerge_LayerCluster_Eta_perlayer" + istr1,
0600                      "Num Merge LayerCluster Eta per Layer Cluster for layer " + istr2,
0601                      nintEta_,
0602                      minEta_,
0603                      maxEta_);
0604     histograms.h_denom_layercl_eta_perlayer[ilayer] =
0605         ibook.book1D("Denom_LayerCluster_Eta_perlayer" + istr1,
0606                      "Denom LayerCluster Eta per Layer Cluster for layer " + istr2,
0607                      nintEta_,
0608                      minEta_,
0609                      maxEta_);
0610     histograms.h_num_layercl_phi_perlayer[ilayer] =
0611         ibook.book1D("Num_LayerCluster_Phi_perlayer" + istr1,
0612                      "Num LayerCluster Phi per Layer Cluster for layer " + istr2,
0613                      nintPhi_,
0614                      minPhi_,
0615                      maxPhi_);
0616     histograms.h_numMerge_layercl_phi_perlayer[ilayer] =
0617         ibook.book1D("NumMerge_LayerCluster_Phi_perlayer" + istr1,
0618                      "Num Merge LayerCluster Phi per Layer Cluster for layer " + istr2,
0619                      nintPhi_,
0620                      minPhi_,
0621                      maxPhi_);
0622     histograms.h_denom_layercl_phi_perlayer[ilayer] =
0623         ibook.book1D("Denom_LayerCluster_Phi_perlayer" + istr1,
0624                      "Denom LayerCluster Phi per Layer Cluster for layer " + istr2,
0625                      nintPhi_,
0626                      minPhi_,
0627                      maxPhi_);
0628   }
0629   //---------------------------------------------------------------------------------------------------------------------------
0630 }
0631 
0632 void BarrelVHistoProducerAlgo::bookClusterHistos_CellLevel(DQMStore::IBooker& ibook,
0633                                                            Histograms& histograms,
0634                                                            unsigned int layers) {
0635   //----------------------------------------------------------------------------------------------------------------------------
0636   for (unsigned ilayer = 0; ilayer < layers; ++ilayer) {
0637     auto istr1 = std::to_string(ilayer);
0638     std::string istr2 = istr1 + " in barrel";
0639     istr1 = "_barrel" + istr1;
0640     histograms.h_cellAssociation_perlayer[ilayer] =
0641         ibook.book1D("cellAssociation_perlayer" + istr1, "Cell Association for layer " + istr2, 5, -4., 1.);
0642     histograms.h_cellAssociation_perlayer[ilayer]->setBinLabel(2, "TN(purity)");
0643     histograms.h_cellAssociation_perlayer[ilayer]->setBinLabel(3, "FN(ineff.)");
0644     histograms.h_cellAssociation_perlayer[ilayer]->setBinLabel(4, "FP(fake)");
0645     histograms.h_cellAssociation_perlayer[ilayer]->setBinLabel(5, "TP(eff.)");
0646   }
0647 }
0648 //----------------------------------------------------------------------------------------------------------------------------
0649 
0650 void BarrelVHistoProducerAlgo::fill_info_histos(const Histograms& histograms, unsigned int layers) const {
0651   // Save some info straight from geometry to avoid mistakes from updates
0652   //----------- TODO ----------------------------------------------------------
0653   // For now values returned for 'lastLayerFHzp': '104', 'lastLayerFHzm': '52' are not the one expected.
0654   // Will come back to this when there will be info in CMSSW to put in DQM file.
0655   histograms.lastLayerEB->Fill(0);
0656   histograms.lastLayerHB->Fill(recHitTools_->lastLayerBarrel());
0657 }
0658 
0659 void BarrelVHistoProducerAlgo::fill_caloparticle_histos(
0660     const Histograms& histograms,
0661     int pdgid,
0662     const CaloParticle& caloParticle,
0663     std::vector<SimVertex> const& simVertices,
0664     unsigned int layers,
0665     std::unordered_map<DetId, const unsigned int> const& barrelHitMap,
0666     MultiVectorManager<reco::PFRecHit> const& barrelHits) const {
0667   const auto eta = getEta(caloParticle.eta());
0668   if (histograms.h_caloparticle_eta.count(pdgid)) {
0669     histograms.h_caloparticle_eta.at(pdgid)->Fill(eta);
0670   }
0671   if (histograms.h_caloparticle_eta_Zorigin.count(pdgid)) {
0672     histograms.h_caloparticle_eta_Zorigin.at(pdgid)->Fill(
0673         simVertices.at(caloParticle.g4Tracks()[0].vertIndex()).position().z(), eta);
0674   }
0675 
0676   if (histograms.h_caloparticle_energy.count(pdgid)) {
0677     histograms.h_caloparticle_energy.at(pdgid)->Fill(caloParticle.energy());
0678   }
0679   if (histograms.h_caloparticle_pt.count(pdgid)) {
0680     histograms.h_caloparticle_pt.at(pdgid)->Fill(caloParticle.pt());
0681   }
0682   if (histograms.h_caloparticle_phi.count(pdgid)) {
0683     histograms.h_caloparticle_phi.at(pdgid)->Fill(caloParticle.phi());
0684   }
0685 
0686   if (histograms.h_caloparticle_nSimClusters.count(pdgid)) {
0687     histograms.h_caloparticle_nSimClusters.at(pdgid)->Fill(caloParticle.simClusters().size());
0688 
0689     int simHits = 0;
0690     int minLayerId = 999;
0691     int maxLayerId = 0;
0692 
0693     int simHits_matched = 0;
0694     int minLayerId_matched = 999;
0695     int maxLayerId_matched = 0;
0696 
0697     float energy = 0.;
0698     std::map<int, double> totenergy_layer;
0699 
0700     float hitEnergyWeight_invSum = 0;
0701     std::vector<std::pair<DetId, float>> haf_cp;
0702     for (const auto& sc : caloParticle.simClusters()) {
0703       LogDebug("BarrelValidator") << " This sim cluster has " << sc->hits_and_fractions().size() << " simHits and "
0704                                   << sc->energy() << " energy. " << std::endl;
0705       simHits += sc->hits_and_fractions().size();
0706       for (auto const& h_and_f : sc->hits_and_fractions()) {
0707         const auto hitDetId = h_and_f.first;
0708         const int layerId = recHitTools_->getLayerWithOffset(hitDetId);
0709         // set to 0 if matched RecHit not found
0710         int layerId_matched_min = 999;
0711         int layerId_matched_max = 0;
0712         std::unordered_map<DetId, const unsigned int>::const_iterator itcheck = barrelHitMap.find(hitDetId);
0713         if (itcheck != barrelHitMap.end()) {
0714           layerId_matched_min = layerId;
0715           layerId_matched_max = layerId;
0716           simHits_matched++;
0717 
0718           const auto hitEn = (barrelHits[itcheck->second]).energy();
0719           hitEnergyWeight_invSum += pow(hitEn, 2);
0720           const auto hitFr = h_and_f.second;
0721           const auto hitEnFr = hitEn * hitFr;
0722           energy += hitEnFr;
0723           histograms.h_caloparticle_nHits_matched_energy.at(pdgid)->Fill(hitEnFr);
0724           histograms.h_caloparticle_nHits_matched_energy_layer.at(pdgid)->Fill(layerId, hitEnFr);
0725 
0726           if (totenergy_layer.find(layerId) != totenergy_layer.end()) {
0727             totenergy_layer[layerId] = totenergy_layer.at(layerId) + hitEn;
0728           } else {
0729             totenergy_layer.emplace(layerId, hitEn);
0730           }
0731           if (caloParticle.simClusters().size() == 1)
0732             histograms.h_caloparticle_nHits_matched_energy_layer_1SimCl.at(pdgid)->Fill(layerId, hitEnFr);
0733 
0734           auto found = std::find_if(std::begin(haf_cp),
0735                                     std::end(haf_cp),
0736                                     [&hitDetId](const std::pair<DetId, float>& v) { return v.first == hitDetId; });
0737           if (found != haf_cp.end())
0738             found->second += hitFr;
0739           else
0740             haf_cp.emplace_back(hitDetId, hitFr);
0741 
0742         } else {
0743           LogDebug("BarrelValidator") << "   matched to RecHit NOT found !" << std::endl;
0744         }
0745 
0746         minLayerId = std::min(minLayerId, layerId);
0747         maxLayerId = std::max(maxLayerId, layerId);
0748         minLayerId_matched = std::min(minLayerId_matched, layerId_matched_min);
0749         maxLayerId_matched = std::max(maxLayerId_matched, layerId_matched_max);
0750       }
0751       LogDebug("BarrelValidator") << std::endl;
0752     }  // End loop over SimClusters of CaloParticle
0753     if (hitEnergyWeight_invSum)
0754       hitEnergyWeight_invSum = 1 / hitEnergyWeight_invSum;
0755 
0756     histograms.h_caloparticle_firstlayer.at(pdgid)->Fill(minLayerId);
0757     histograms.h_caloparticle_lastlayer.at(pdgid)->Fill(maxLayerId);
0758     histograms.h_caloparticle_layersnum.at(pdgid)->Fill(int(maxLayerId - minLayerId));
0759 
0760     histograms.h_caloparticle_firstlayer_matchedtoRecHit.at(pdgid)->Fill(minLayerId_matched);
0761     histograms.h_caloparticle_lastlayer_matchedtoRecHit.at(pdgid)->Fill(maxLayerId_matched);
0762     histograms.h_caloparticle_layersnum_matchedtoRecHit.at(pdgid)->Fill(int(maxLayerId_matched - minLayerId_matched));
0763 
0764     histograms.h_caloparticle_nHitsInSimClusters.at(pdgid)->Fill((float)simHits);
0765     histograms.h_caloparticle_nHitsInSimClusters_matchedtoRecHit.at(pdgid)->Fill((float)simHits_matched);
0766     histograms.h_caloparticle_selfenergy.at(pdgid)->Fill((float)energy);
0767     histograms.h_caloparticle_energyDifference.at(pdgid)->Fill((float)1. - energy / caloParticle.energy());
0768 
0769     //Calculate sum energy per-layer
0770     auto i = totenergy_layer.begin();
0771     double sum_energy = 0.0;
0772     while (i != totenergy_layer.end()) {
0773       sum_energy += i->second;
0774       histograms.h_caloparticle_sum_energy_layer.at(pdgid)->Fill(i->first, sum_energy / caloParticle.energy() * 100.);
0775       i++;
0776     }
0777 
0778     for (auto const& haf : haf_cp) {
0779       if (!recHitTools_->isBarrel(haf.first))
0780         continue;
0781       const auto hitEn = (barrelHits[barrelHitMap.find(haf.first)->second]).energy();
0782       const auto weight = pow(hitEn, 2);
0783       histograms.h_caloparticle_fractions.at(pdgid)->Fill(haf.second, weight * hitEnergyWeight_invSum);
0784       histograms.h_caloparticle_fractions_weight.at(pdgid)->Fill(haf.second, weight * hitEnergyWeight_invSum, weight);
0785     }
0786   }
0787 }
0788 
0789 void BarrelVHistoProducerAlgo::BarrelVHistoProducerAlgo::fill_simCluster_histos(
0790     const Histograms& histograms, std::vector<SimCluster> const& simClusters, unsigned int layers) const {
0791   //To keep track of total num of simClusters per layer
0792   //tnscpl[layerid]
0793   std::vector<int> tnscpl(1000, 0);  //tnscpl.clear(); tnscpl.reserve(1000);
0794 
0795   //loop through simClusters
0796   for (const auto& sc : simClusters) {
0797     //To keep track if we added the simCluster in a specific layer
0798     std::vector<int> occurenceSCinlayer(1000, 0);  //[layerid][0 if not added]
0799 
0800     //loop through hits of the simCluster
0801     for (const auto& hAndF : sc.hits_and_fractions()) {
0802       const DetId sh_detid = hAndF.first;
0803 
0804       if (recHitTools_->isBarrel(sh_detid)) {
0805         //The layer the cluster belongs to. As mentioned in the mapping above, it takes into account -z and +z.
0806         int layerid = recHitTools_->getLayerWithOffset(sh_detid);
0807         //zside that the current cluster belongs to.
0808         if (occurenceSCinlayer[layerid] == 0) {
0809           tnscpl[layerid]++;
0810         }
0811         occurenceSCinlayer[layerid]++;
0812       }
0813     }  //end of loop through hits
0814   }  //end of loop through SimClusters of the event
0815 
0816   //Per layer : Loop 0->99
0817   for (unsigned ilayer = 0; ilayer < layers; ++ilayer) {
0818     if (histograms.h_simclusternum_perlayer.count(ilayer))
0819       histograms.h_simclusternum_perlayer.at(ilayer)->Fill(tnscpl[ilayer]);
0820   }
0821 }
0822 
0823 void BarrelVHistoProducerAlgo::BarrelVHistoProducerAlgo::fill_simClusterAssociation_histos(
0824     const Histograms& histograms,
0825     const int count,
0826     edm::Handle<reco::CaloClusterCollection> clusterHandle,
0827     const reco::CaloClusterCollection& clusters,
0828     edm::Handle<std::vector<SimCluster>> simClusterHandle,
0829     std::vector<SimCluster> const& simClusters,
0830     std::vector<size_t> const& sCIndices,
0831     const std::vector<float>& mask,
0832     std::unordered_map<DetId, const unsigned int> const& barrelHitMap,
0833     unsigned int layers,
0834     const ticl::RecoToSimCollectionWithSimClusters& scsInLayerClusterMap,
0835     const ticl::SimToRecoCollectionWithSimClusters& lcsInSimClusterMap,
0836     MultiVectorManager<reco::PFRecHit> const& barrelHits) const {
0837   //Each event to be treated as two events: an event in +ve endcap,
0838   //plus another event in -ve endcap. In this spirit there will be
0839   //a layer variable (layerid) that maps the layers in :
0840   //-z: 0->49
0841   //+z: 50->99
0842 
0843   //Will add some general plots on the specific mask in the future.
0844 
0845   layerClusters_to_SimClusters(histograms,
0846                                count,
0847                                clusterHandle,
0848                                clusters,
0849                                simClusterHandle,
0850                                simClusters,
0851                                sCIndices,
0852                                mask,
0853                                barrelHitMap,
0854                                layers,
0855                                scsInLayerClusterMap,
0856                                lcsInSimClusterMap,
0857                                barrelHits);
0858 }
0859 
0860 void BarrelVHistoProducerAlgo::fill_cluster_histos(const Histograms& histograms,
0861                                                    const int count,
0862                                                    const reco::CaloCluster& cluster) const {
0863   if (!recHitTools_->isBarrel(cluster.seed()))
0864     return;
0865   const auto eta = getEta(cluster.eta());
0866   histograms.h_cluster_eta[count]->Fill(eta);
0867 }
0868 
0869 void BarrelVHistoProducerAlgo::layerClusters_to_CaloParticles(
0870     const Histograms& histograms,
0871     edm::Handle<reco::CaloClusterCollection> clusterHandle,
0872     const reco::CaloClusterCollection& clusters,
0873     edm::Handle<std::vector<CaloParticle>> caloParticleHandle,
0874     std::vector<CaloParticle> const& cP,
0875     std::vector<size_t> const& cPIndices,
0876     std::vector<size_t> const& cPSelectedIndices,
0877     std::unordered_map<DetId, const unsigned int> const& barrelHitMap,
0878     unsigned int layers,
0879     const ticl::RecoToSimCollection& cpsInLayerClusterMap,
0880     const ticl::SimToRecoCollection& cPOnLayerMap,
0881     MultiVectorManager<reco::PFRecHit> const& barrelHits) const {
0882   const auto nLayerClusters = clusters.size();
0883 
0884   std::unordered_map<DetId, std::vector<BarrelVHistoProducerAlgo::detIdInfoInCluster>> detIdToCaloParticleId_Map;
0885   std::unordered_map<DetId, std::vector<BarrelVHistoProducerAlgo::detIdInfoInCluster>> detIdToLayerClusterId_Map;
0886 
0887   // The association has to be done in an all-vs-all fashion.
0888   // For this reason use the full set of CaloParticles, with the only filter on bx
0889   for (const auto& cpId : cPIndices) {
0890     for (const auto& simCluster : cP[cpId].simClusters()) {
0891       for (const auto& it_haf : simCluster->hits_and_fractions()) {
0892         const DetId hitid = (it_haf.first);
0893         if (barrelHitMap.find(hitid) != barrelHitMap.end()) {
0894           if (detIdToCaloParticleId_Map.find(hitid) == detIdToCaloParticleId_Map.end()) {
0895             detIdToCaloParticleId_Map[hitid] = std::vector<BarrelVHistoProducerAlgo::detIdInfoInCluster>();
0896             detIdToCaloParticleId_Map[hitid].emplace_back(
0897                 BarrelVHistoProducerAlgo::detIdInfoInCluster{cpId, it_haf.second});
0898           } else {
0899             auto findHitIt =
0900                 std::find(detIdToCaloParticleId_Map[hitid].begin(),
0901                           detIdToCaloParticleId_Map[hitid].end(),
0902                           BarrelVHistoProducerAlgo::detIdInfoInCluster{
0903                               cpId, 0.f});  // only the first element is used for the matching (overloaded operator==)
0904             if (findHitIt != detIdToCaloParticleId_Map[hitid].end())
0905               findHitIt->fraction += it_haf.second;
0906             else
0907               detIdToCaloParticleId_Map[hitid].emplace_back(
0908                   BarrelVHistoProducerAlgo::detIdInfoInCluster{cpId, it_haf.second});
0909           }
0910         }
0911       }
0912     }
0913   }
0914 
0915   for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) {
0916     const auto& hits_and_fractions = clusters[lcId].hitsAndFractions();
0917     const auto numberOfHitsInLC = hits_and_fractions.size();
0918 
0919     // This vector will store, for each hit in the Layercluster, the index of
0920     // the CaloParticle that contributed the most, in terms of energy, to it.
0921     // Special values are:
0922     //
0923     // -2  --> the reconstruction fraction of the RecHit is 0 (used in the past to monitor Halo Hits)
0924     // -3  --> same as before with the added condition that no CaloParticle has been linked to this RecHit
0925     // -1  --> the reco fraction is >0, but no CaloParticle has been linked to it
0926     // >=0 --> index of the linked CaloParticle
0927     std::vector<int> hitsToCaloParticleId(numberOfHitsInLC);
0928     const auto firstHitDetId = hits_and_fractions[0].first;
0929     int lcLayerId = recHitTools_->getLayerWithOffset(firstHitDetId);
0930     bool isBarrel = recHitTools_->isBarrel(firstHitDetId);
0931     if (!isBarrel)
0932       continue;
0933     // This will store the fraction of the CaloParticle energy shared with the LayerCluster: e_shared/cp_energy
0934     std::unordered_map<unsigned, float> CPEnergyInLC;
0935 
0936     for (unsigned int iHit = 0; iHit < numberOfHitsInLC; iHit++) {
0937       const DetId rh_detid = hits_and_fractions[iHit].first;
0938       const auto rhFraction = hits_and_fractions[iHit].second;
0939 
0940       float hit_energy;
0941       if (!isBarrel) {
0942         continue;
0943       } else {
0944         std::unordered_map<DetId, const unsigned int>::const_iterator itcheck = barrelHitMap.find(rh_detid);
0945         const reco::PFRecHit* hit = &(barrelHits[itcheck->second]);
0946         hit_energy = hit->energy();
0947       }
0948       if (detIdToLayerClusterId_Map.find(rh_detid) == detIdToLayerClusterId_Map.end()) {
0949         detIdToLayerClusterId_Map[rh_detid] = std::vector<BarrelVHistoProducerAlgo::detIdInfoInCluster>();
0950       }
0951       detIdToLayerClusterId_Map[rh_detid].emplace_back(BarrelVHistoProducerAlgo::detIdInfoInCluster{lcId, rhFraction});
0952 
0953       const auto& hit_find_in_CP = detIdToCaloParticleId_Map.find(rh_detid);
0954 
0955       // if the fraction is zero or the hit does not belong to any calo
0956       // particle, set the caloparticleId for the hit to -1 this will
0957       // contribute to the number of noise hits
0958 
0959       // MR Remove the case in which the fraction is 0, since this could be a
0960       // real hit that has been marked as halo.
0961       if (rhFraction == 0.) {
0962         hitsToCaloParticleId[iHit] = -2;
0963       }
0964       if (hit_find_in_CP == detIdToCaloParticleId_Map.end()) {
0965         hitsToCaloParticleId[iHit] -= 1;
0966       } else {
0967         auto maxCPEnergyInLC = 0.f;
0968         auto maxCPId = -1;
0969         for (auto& h : hit_find_in_CP->second) {
0970           const auto iCP = h.clusterId;
0971           CPEnergyInLC[iCP] += h.fraction * hit_energy;
0972           // Keep track of which CaloParticle contributed the most, in terms
0973           // of energy, to this specific LayerCluster.
0974           if (CPEnergyInLC[iCP] > maxCPEnergyInLC) {
0975             maxCPEnergyInLC = CPEnergyInLC[iCP];
0976             maxCPId = iCP;
0977           }
0978         }
0979         hitsToCaloParticleId[iHit] = maxCPId;
0980       }
0981       histograms.h_cellAssociation_perlayer.at(lcLayerId)->Fill(
0982           hitsToCaloParticleId[iHit] > 0. ? 0. : hitsToCaloParticleId[iHit]);
0983     }  // End loop over hits on a LayerCluster
0984 
0985   }  // End of loop over LayerClusters
0986 
0987   // Fill the plots to compute the different metrics linked to
0988   // reco-level, namely fake-rate an merge-rate. In this loop should *not*
0989   // restrict only to the selected caloParticles.
0990   for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) {
0991     const auto firstHitDetId = (clusters[lcId].hitsAndFractions())[0].first;
0992     bool isBarrel = recHitTools_->isBarrel(firstHitDetId);
0993     int lcLayerId = recHitTools_->getLayerWithOffset(firstHitDetId);
0994     if (!isBarrel)
0995       continue;
0996     histograms.h_denom_layercl_eta_perlayer.at(lcLayerId)->Fill(clusters[lcId].eta());
0997     histograms.h_denom_layercl_phi_perlayer.at(lcLayerId)->Fill(clusters[lcId].phi());
0998     //
0999     const edm::Ref<reco::CaloClusterCollection> lcRef(clusterHandle, lcId);
1000     const auto& cpsIt = cpsInLayerClusterMap.find(lcRef);
1001     if (cpsIt == cpsInLayerClusterMap.end())
1002       continue;
1003 
1004     const auto lc_en = clusters[lcId].energy();
1005     const auto& cps = cpsIt->val;
1006     if (lc_en == 0. && !cps.empty()) {
1007       for (const auto& cpPair : cps)
1008         histograms.h_score_layercl2caloparticle_perlayer.at(lcLayerId)->Fill(cpPair.second);
1009       continue;
1010     }
1011     for (const auto& cpPair : cps) {
1012       LogDebug("BarrelValidator") << "layerCluster Id: \t" << lcId << "\t CP id: \t" << cpPair.first.index()
1013                                   << "\t score \t" << cpPair.second << std::endl;
1014       histograms.h_score_layercl2caloparticle_perlayer.at(lcLayerId)->Fill(cpPair.second);
1015       auto const& cp_linked =
1016           std::find_if(std::begin(cPOnLayerMap[cpPair.first]),
1017                        std::end(cPOnLayerMap[cpPair.first]),
1018                        [&lcRef](const std::pair<edm::Ref<reco::CaloClusterCollection>, std::pair<float, float>>& p) {
1019                          return p.first == lcRef;
1020                        });
1021       if (cp_linked ==
1022           cPOnLayerMap[cpPair.first].end())  // This should never happen by construction of the association maps
1023         continue;
1024       histograms.h_sharedenergy_layercl2caloparticle_perlayer.at(lcLayerId)->Fill(cp_linked->second.first / lc_en,
1025                                                                                   lc_en);
1026       histograms.h_energy_vs_score_layercl2caloparticle_perlayer.at(lcLayerId)->Fill(cpPair.second,
1027                                                                                      cp_linked->second.first / lc_en);
1028     }
1029     const auto assoc =
1030         std::count_if(std::begin(cps), std::end(cps), [](const auto& obj) { return obj.second < ScoreCutLCtoCP_; });
1031     if (assoc) {
1032       histograms.h_num_layercl_eta_perlayer.at(lcLayerId)->Fill(clusters[lcId].eta());
1033       histograms.h_num_layercl_phi_perlayer.at(lcLayerId)->Fill(clusters[lcId].phi());
1034       if (assoc > 1) {
1035         histograms.h_numMerge_layercl_eta_perlayer.at(lcLayerId)->Fill(clusters[lcId].eta());
1036         histograms.h_numMerge_layercl_phi_perlayer.at(lcLayerId)->Fill(clusters[lcId].phi());
1037       }
1038       const auto& best = std::min_element(
1039           std::begin(cps), std::end(cps), [](const auto& obj1, const auto& obj2) { return obj1.second < obj2.second; });
1040       const auto& best_cp_linked =
1041           std::find_if(std::begin(cPOnLayerMap[best->first]),
1042                        std::end(cPOnLayerMap[best->first]),
1043                        [&lcRef](const std::pair<edm::Ref<reco::CaloClusterCollection>, std::pair<float, float>>& p) {
1044                          return p.first == lcRef;
1045                        });
1046       if (best_cp_linked ==
1047           cPOnLayerMap[best->first].end())  // This should never happen by construction of the association maps
1048         continue;
1049       histograms.h_sharedenergy_layercl2caloparticle_vs_eta_perlayer.at(lcLayerId)->Fill(
1050           clusters[lcId].eta(), best_cp_linked->second.first / lc_en);
1051       histograms.h_sharedenergy_layercl2caloparticle_vs_phi_perlayer.at(lcLayerId)->Fill(
1052           clusters[lcId].phi(), best_cp_linked->second.first / lc_en);
1053     }
1054   }  // End of loop over LayerClusters
1055 
1056   // Here Fill the plots to compute the different metrics linked to
1057   // gen-level, namely efficiency and duplicate. In this loop should restrict
1058   // only to the selected caloParaticles.
1059   for (const auto& cpId : cPSelectedIndices) {
1060     const edm::Ref<CaloParticleCollection> cpRef(caloParticleHandle, cpId);
1061     const auto& lcsIt = cPOnLayerMap.find(cpRef);
1062 
1063     std::map<unsigned int, float> cPEnergyOnLayer;
1064     for (unsigned int layerId = 0; layerId < layers; ++layerId)
1065       cPEnergyOnLayer[layerId] = 0;
1066 
1067     for (const auto& simCluster : cP[cpId].simClusters()) {
1068       for (const auto& it_haf : simCluster->hits_and_fractions()) {
1069         const DetId hitid = (it_haf.first);
1070         bool isBarrel = recHitTools_->isBarrel(hitid);
1071         auto hitLayerId = recHitTools_->getLayerWithOffset(hitid);
1072         if (!isBarrel) {
1073           continue;
1074         } else {
1075           std::unordered_map<DetId, const unsigned int>::const_iterator itcheck = barrelHitMap.find(hitid);
1076           if (itcheck != barrelHitMap.end()) {
1077             const reco::PFRecHit* hit = &(barrelHits[itcheck->second]);
1078             cPEnergyOnLayer[hitLayerId] += it_haf.second * hit->energy();
1079           }
1080         }
1081       }
1082     }
1083 
1084     for (unsigned int layerId = 0; layerId < layers; ++layerId) {
1085       if (!cPEnergyOnLayer[layerId])
1086         continue;
1087 
1088       histograms.h_denom_caloparticle_eta_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().eta());
1089       histograms.h_denom_caloparticle_phi_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().phi());
1090 
1091       if (lcsIt == cPOnLayerMap.end())
1092         continue;
1093       const auto& lcs = lcsIt->val;
1094 
1095       auto getLCLayerId = [&](const unsigned int lcId) {
1096         const auto firstHitDetId = (clusters[lcId].hitsAndFractions())[0].first;
1097         bool isBarrel = recHitTools_->isBarrel(firstHitDetId);
1098         if (!isBarrel)
1099           return 9999u;
1100         auto lcLayerId = recHitTools_->getLayerWithOffset(firstHitDetId);
1101         return lcLayerId;
1102       };
1103 
1104       for (const auto& lcPair : lcs) {
1105         if (getLCLayerId(lcPair.first.index()) != layerId)
1106           continue;
1107         histograms.h_score_caloparticle2layercl_perlayer.at(layerId)->Fill(lcPair.second.second);
1108         histograms.h_sharedenergy_caloparticle2layercl_perlayer.at(layerId)->Fill(
1109             lcPair.second.first / cPEnergyOnLayer[layerId], cPEnergyOnLayer[layerId]);
1110         histograms.h_energy_vs_score_caloparticle2layercl_perlayer.at(layerId)->Fill(
1111             lcPair.second.second, lcPair.second.first / cPEnergyOnLayer[layerId]);
1112       }
1113       const auto assoc = std::count_if(std::begin(lcs), std::end(lcs), [&](const auto& obj) {
1114         if (getLCLayerId(obj.first.index()) != layerId)
1115           return false;
1116         else
1117           return obj.second.second < ScoreCutCPtoLC_;
1118       });
1119       if (assoc) {
1120         histograms.h_num_caloparticle_eta_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().eta());
1121         histograms.h_num_caloparticle_phi_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().phi());
1122         if (assoc > 1) {
1123           histograms.h_numDup_caloparticle_eta_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().eta());
1124           histograms.h_numDup_caloparticle_phi_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().phi());
1125         }
1126         const auto best = std::min_element(std::begin(lcs), std::end(lcs), [&](const auto& obj1, const auto& obj2) {
1127           if (getLCLayerId(obj1.first.index()) != layerId)
1128             return false;
1129           else if (getLCLayerId(obj2.first.index()) == layerId)
1130             return obj1.second.second < obj2.second.second;
1131           else
1132             return true;
1133         });
1134         histograms.h_sharedenergy_caloparticle2layercl_vs_eta_perlayer.at(layerId)->Fill(
1135             cP[cpId].g4Tracks()[0].momentum().eta(), best->second.first / cPEnergyOnLayer[layerId]);
1136         histograms.h_sharedenergy_caloparticle2layercl_vs_phi_perlayer.at(layerId)->Fill(
1137             cP[cpId].g4Tracks()[0].momentum().phi(), best->second.first / cPEnergyOnLayer[layerId]);
1138       }
1139     }
1140   }
1141 }
1142 
1143 void BarrelVHistoProducerAlgo::layerClusters_to_SimClusters(
1144     const Histograms& histograms,
1145     const int count,
1146     edm::Handle<reco::CaloClusterCollection> clusterHandle,
1147     const reco::CaloClusterCollection& clusters,
1148     edm::Handle<std::vector<SimCluster>> simClusterHandle,
1149     std::vector<SimCluster> const& sC,
1150     std::vector<size_t> const& sCIndices,
1151     const std::vector<float>& mask,
1152     std::unordered_map<DetId, const unsigned int> const& barrelHitMap,
1153     unsigned int layers,
1154     const ticl::RecoToSimCollectionWithSimClusters& scsInLayerClusterMap,
1155     const ticl::SimToRecoCollectionWithSimClusters& lcsInSimClusterMap,
1156     MultiVectorManager<reco::PFRecHit> const& barrelHits) const {
1157   // Here fill the plots to compute the different metrics linked to
1158   // reco-level, namely fake-rate and merge-rate. In this loop should *not*
1159   // restrict only to the selected SimClusters.
1160   for (unsigned int lcId = 0; lcId < clusters.size(); ++lcId) {
1161     if (mask[lcId] != 0.) {
1162       LogDebug("BarrelValidator") << "Skipping layer cluster " << lcId << " not belonging to mask" << std::endl;
1163       continue;
1164     }
1165     const auto firstHitDetId = (clusters[lcId].hitsAndFractions())[0].first;
1166     bool isBarrel = recHitTools_->isBarrel(firstHitDetId);
1167     auto lcLayerId = recHitTools_->getLayerWithOffset(firstHitDetId);
1168     if (!isBarrel)
1169       continue;
1170     //Although the ones below are already created in the LC to CP association, will
1171     //recreate them here since in the post processor it looks in a specific directory.
1172     histograms.h_denom_layercl_in_simcl_eta_perlayer[count].at(lcLayerId)->Fill(clusters[lcId].eta());
1173     histograms.h_denom_layercl_in_simcl_phi_perlayer[count].at(lcLayerId)->Fill(clusters[lcId].phi());
1174     //
1175     const edm::Ref<reco::CaloClusterCollection> lcRef(clusterHandle, lcId);
1176     const auto& scsIt = scsInLayerClusterMap.find(lcRef);
1177     if (scsIt == scsInLayerClusterMap.end())
1178       continue;
1179 
1180     const auto lc_en = clusters[lcId].energy();
1181     const auto& scs = scsIt->val;
1182     // If a reconstructed LayerCluster has energy 0 but is linked to at least a
1183     // SimCluster, then his score should be 1 as set in the associator
1184     if (lc_en == 0. && !scs.empty()) {
1185       for (const auto& scPair : scs) {
1186         histograms.h_score_layercl2simcluster_perlayer[count].at(lcLayerId)->Fill(scPair.second);
1187       }
1188       continue;
1189     }
1190     //Loop through all SimClusters linked to the layer cluster under study
1191     for (const auto& scPair : scs) {
1192       LogDebug("BarrelValidator") << "layerCluster Id: \t" << lcId << "\t SC id: \t" << scPair.first.index()
1193                                   << "\t score \t" << scPair.second << std::endl;
1194       //This should be filled #layerClusters in layer x #linked SimClusters
1195       histograms.h_score_layercl2simcluster_perlayer[count].at(lcLayerId)->Fill(scPair.second);
1196       auto const& sc_linked =
1197           std::find_if(std::begin(lcsInSimClusterMap[scPair.first]),
1198                        std::end(lcsInSimClusterMap[scPair.first]),
1199                        [&lcRef](const std::pair<edm::Ref<reco::CaloClusterCollection>, std::pair<float, float>>& p) {
1200                          return p.first == lcRef;
1201                        });
1202       if (sc_linked ==
1203           lcsInSimClusterMap[scPair.first].end())  // This should never happen by construction of the association maps
1204         continue;
1205       histograms.h_sharedenergy_layercl2simcluster_perlayer[count].at(lcLayerId)->Fill(sc_linked->second.first / lc_en,
1206                                                                                        lc_en);
1207       histograms.h_energy_vs_score_layercl2simcluster_perlayer[count].at(lcLayerId)->Fill(
1208           scPair.second, sc_linked->second.first / lc_en);
1209     }
1210     //Here he counts how many of the linked SimClusters of the layer cluster under study have a score above a certain value.
1211     const auto assoc =
1212         std::count_if(std::begin(scs), std::end(scs), [](const auto& obj) { return obj.second < ScoreCutLCtoSC_; });
1213     if (assoc) {
1214       histograms.h_num_layercl_in_simcl_eta_perlayer[count].at(lcLayerId)->Fill(clusters[lcId].eta());
1215       histograms.h_num_layercl_in_simcl_phi_perlayer[count].at(lcLayerId)->Fill(clusters[lcId].phi());
1216       if (assoc > 1) {
1217         histograms.h_numMerge_layercl_in_simcl_eta_perlayer[count].at(lcLayerId)->Fill(clusters[lcId].eta());
1218         histograms.h_numMerge_layercl_in_simcl_phi_perlayer[count].at(lcLayerId)->Fill(clusters[lcId].phi());
1219       }
1220       const auto& best = std::min_element(
1221           std::begin(scs), std::end(scs), [](const auto& obj1, const auto& obj2) { return obj1.second < obj2.second; });
1222       //From all SimClusters he founds the one with the best (lowest) score and takes his scId
1223       const auto& best_sc_linked =
1224           std::find_if(std::begin(lcsInSimClusterMap[best->first]),
1225                        std::end(lcsInSimClusterMap[best->first]),
1226                        [&lcRef](const std::pair<edm::Ref<reco::CaloClusterCollection>, std::pair<float, float>>& p) {
1227                          return p.first == lcRef;
1228                        });
1229       if (best_sc_linked ==
1230           lcsInSimClusterMap[best->first].end())  // This should never happen by construction of the association maps
1231         continue;
1232       histograms.h_sharedenergy_layercl2simcluster_vs_eta_perlayer[count].at(lcLayerId)->Fill(
1233           clusters[lcId].eta(), best_sc_linked->second.first / lc_en);
1234       histograms.h_sharedenergy_layercl2simcluster_vs_phi_perlayer[count].at(lcLayerId)->Fill(
1235           clusters[lcId].phi(), best_sc_linked->second.first / lc_en);
1236     }
1237   }  // End of loop over LayerClusters
1238 
1239   // Fill the plots to compute the different metrics linked to
1240   // gen-level, namely efficiency and duplicate. In this loop should restrict
1241   // only to the selected SimClusters.
1242   for (const auto& scId : sCIndices) {
1243     const edm::Ref<SimClusterCollection> scRef(simClusterHandle, scId);
1244     const auto& lcsIt = lcsInSimClusterMap.find(scRef);
1245 
1246     std::map<unsigned int, float> sCEnergyOnLayer;
1247     for (unsigned int layerId = 0; layerId < layers * 2; ++layerId)
1248       sCEnergyOnLayer[layerId] = 0;
1249 
1250     for (const auto& it_haf : sC[scId].hits_and_fractions()) {
1251       const DetId hitid = (it_haf.first);
1252       bool isBarrel = recHitTools_->isBarrel(hitid);
1253       auto scLayerId = recHitTools_->getLayerWithOffset(hitid);
1254       if (!isBarrel) {
1255         continue;
1256       } else {
1257         std::unordered_map<DetId, const unsigned int>::const_iterator itcheck = barrelHitMap.find(hitid);
1258         if (itcheck != barrelHitMap.end()) {
1259           const reco::PFRecHit* hit = &(barrelHits[itcheck->second]);
1260           sCEnergyOnLayer[scLayerId] += it_haf.second * hit->energy();
1261         }
1262       }
1263     }
1264 
1265     for (unsigned int layerId = 0; layerId < layers * 2; ++layerId) {
1266       if (!sCEnergyOnLayer[layerId])
1267         continue;
1268 
1269       histograms.h_denom_simcluster_eta_perlayer[count].at(layerId)->Fill(sC[scId].eta());
1270       histograms.h_denom_simcluster_phi_perlayer[count].at(layerId)->Fill(sC[scId].phi());
1271 
1272       if (lcsIt == lcsInSimClusterMap.end())
1273         continue;
1274       const auto& lcs = lcsIt->val;
1275 
1276       auto getLCLayerId = [&](const unsigned int lcId) {
1277         const auto firstHitDetId = (clusters[lcId].hitsAndFractions())[0].first;
1278         bool isBarrel = recHitTools_->isBarrel(firstHitDetId);
1279         if (!isBarrel)
1280           return 9999u;
1281         unsigned int lcLayerId = recHitTools_->getLayerWithOffset(firstHitDetId);
1282         return lcLayerId;
1283       };
1284 
1285       //Loop through layer clusters linked to the SimCluster under study
1286       for (const auto& lcPair : lcs) {
1287         auto lcId = lcPair.first.index();
1288         if (mask[lcId] != 0.) {
1289           LogDebug("BarrelValidator") << "Skipping layer cluster " << lcId << " not belonging to mask" << std::endl;
1290           continue;
1291         }
1292 
1293         if (getLCLayerId(lcId) != layerId)
1294           continue;
1295         histograms.h_score_simcluster2layercl_perlayer[count].at(layerId)->Fill(lcPair.second.second);
1296         histograms.h_sharedenergy_simcluster2layercl_perlayer[count].at(layerId)->Fill(
1297             lcPair.second.first / sCEnergyOnLayer[layerId], sCEnergyOnLayer[layerId]);
1298         histograms.h_energy_vs_score_simcluster2layercl_perlayer[count].at(layerId)->Fill(
1299             lcPair.second.second, lcPair.second.first / sCEnergyOnLayer[layerId]);
1300       }
1301       const auto assoc = std::count_if(std::begin(lcs), std::end(lcs), [&](const auto& obj) {
1302         if (getLCLayerId(obj.first.index()) != layerId)
1303           return false;
1304         else
1305           return obj.second.second < ScoreCutSCtoLC_;
1306       });
1307       if (assoc) {
1308         histograms.h_num_simcluster_eta_perlayer[count].at(layerId)->Fill(sC[scId].eta());
1309         histograms.h_num_simcluster_phi_perlayer[count].at(layerId)->Fill(sC[scId].phi());
1310         if (assoc > 1) {
1311           histograms.h_numDup_simcluster_eta_perlayer[count].at(layerId)->Fill(sC[scId].eta());
1312           histograms.h_numDup_simcluster_phi_perlayer[count].at(layerId)->Fill(sC[scId].phi());
1313         }
1314         const auto best = std::min_element(std::begin(lcs), std::end(lcs), [&](const auto& obj1, const auto& obj2) {
1315           if (getLCLayerId(obj1.first.index()) != layerId)
1316             return false;
1317           else if (getLCLayerId(obj2.first.index()) == layerId)
1318             return obj1.second.second < obj2.second.second;
1319           else
1320             return true;
1321         });
1322         histograms.h_sharedenergy_simcluster2layercl_vs_eta_perlayer[count].at(layerId)->Fill(
1323             sC[scId].eta(), best->second.first / sCEnergyOnLayer[layerId]);
1324         histograms.h_sharedenergy_simcluster2layercl_vs_phi_perlayer[count].at(layerId)->Fill(
1325             sC[scId].phi(), best->second.first / sCEnergyOnLayer[layerId]);
1326       }
1327     }
1328   }
1329 }
1330 
1331 void BarrelVHistoProducerAlgo::fill_generic_cluster_histos(
1332     const Histograms& histograms,
1333     const int count,
1334     edm::Handle<reco::CaloClusterCollection> clusterHandle,
1335     const reco::CaloClusterCollection& clusters,
1336     edm::Handle<std::vector<CaloParticle>> caloParticleHandle,
1337     std::vector<CaloParticle> const& cP,
1338     std::vector<size_t> const& cPIndices,
1339     std::vector<size_t> const& cPSelectedIndices,
1340     std::unordered_map<DetId, const unsigned int> const& barrelHitMap,
1341     unsigned int layers,
1342     const ticl::RecoToSimCollection& cpsInLayerClusterMap,
1343     const ticl::SimToRecoCollection& cPOnLayerMap,
1344     MultiVectorManager<reco::PFRecHit> const& barrelHits) const {
1345   //To keep track of total num of layer clusters per layer
1346   //tnlcpl[layerid]
1347   std::vector<int> tnlcpl(1000, 0);  //tnlcpl.clear(); tnlcpl.reserve(1000);
1348 
1349   layerClusters_to_CaloParticles(histograms,
1350                                  clusterHandle,
1351                                  clusters,
1352                                  caloParticleHandle,
1353                                  cP,
1354                                  cPIndices,
1355                                  cPSelectedIndices,
1356                                  barrelHitMap,
1357                                  layers,
1358                                  cpsInLayerClusterMap,
1359                                  cPOnLayerMap,
1360                                  barrelHits);
1361 
1362   std::vector<double> tecpl(1000, 0.0);  //tecpl.clear(); tecpl.reserve(1000);
1363 
1364   // loop through clusters of the event
1365   for (const auto& lcId : clusters) {
1366     const auto seedid = lcId.seed();
1367     if (!recHitTools_->isBarrel(seedid)) {
1368       continue;
1369     }
1370     const float lc_en = lcId.energy();
1371     int layerid = recHitTools_->getLayerWithOffset(seedid);
1372     //Energy clustered per layer
1373     tecpl[layerid] = tecpl[layerid] + lc_en;
1374     tnlcpl[layerid] = tnlcpl[layerid] + 1;
1375 
1376   }  //end of loop through clusters of the event
1377 
1378   for (unsigned ilayer = 0; ilayer < layers; ++ilayer) {
1379     if (histograms.h_clusternum_perlayer.count(ilayer)) {
1380       histograms.h_clusternum_perlayer.at(ilayer)->Fill(tnlcpl[ilayer]);
1381       histograms.h_energyclustered_perlayer.at(ilayer)->Fill(tecpl[ilayer]);
1382     }
1383   }  //end of loop over layers
1384 }
1385 
1386 double BarrelVHistoProducerAlgo::distance2(const double x1,
1387                                            const double y1,
1388                                            const double x2,
1389                                            const double y2) const {  //distance squared
1390   const double dx = x1 - x2;
1391   const double dy = y1 - y2;
1392   return (dx * dx + dy * dy);
1393 }  //distance squaredq
1394 double BarrelVHistoProducerAlgo::distance(const double x1,
1395                                           const double y1,
1396                                           const double x2,
1397                                           const double y2) const {  //2-d distance on the layer (x-y)
1398   return std::sqrt(distance2(x1, y1, x2, y2));
1399 }
1400 
1401 void BarrelVHistoProducerAlgo::setRecHitTools(std::shared_ptr<hgcal::RecHitTools> recHitTools) {
1402   recHitTools_ = recHitTools;
1403 }
1404 
1405 DetId BarrelVHistoProducerAlgo::findmaxhit(const reco::CaloCluster& cluster,
1406                                            std::unordered_map<DetId, const unsigned int> const& hitMap,
1407                                            MultiVectorManager<reco::PFRecHit> const& hits) const {
1408   const auto& hits_and_fractions = cluster.hitsAndFractions();
1409 
1410   DetId themaxid;
1411   double maxene = 0.;
1412   for (std::vector<std::pair<DetId, float>>::const_iterator it_haf = hits_and_fractions.begin();
1413        it_haf != hits_and_fractions.end();
1414        ++it_haf) {
1415     const DetId rh_detid = it_haf->first;
1416     if (!recHitTools_->isBarrel(rh_detid))
1417       continue;
1418     const auto hitEn = (hits[hitMap.find(rh_detid)->second]).energy();
1419     if (maxene < hitEn) {
1420       maxene = hitEn;
1421       themaxid = rh_detid;
1422     }
1423   }
1424 
1425   return themaxid;
1426 }
1427 
1428 double BarrelVHistoProducerAlgo::getEta(double eta) const {
1429   if (useFabsEta_)
1430     return fabs(eta);
1431   else
1432     return eta;
1433 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!