Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-10-10 04:02:18

0001 #include <numeric>
0002 #include <iomanip>
0003 #include <sstream>
0004 
0005 #include "Validation/HGCalValidation/interface/HGVHistoProducerAlgo.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 HGCAL 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 const double ScoreCutTStoSTSFakeMerge_[] = {0.6, FLT_MIN};  //1.e-09
0021 const double ScoreCutSTStoTSPurDup_[] = {0.2, FLT_MIN};     //1.e-11
0022 
0023 HGVHistoProducerAlgo::HGVHistoProducerAlgo(const edm::ParameterSet& pset)
0024     :  //parameters for eta
0025       minEta_(pset.getParameter<double>("minEta")),
0026       maxEta_(pset.getParameter<double>("maxEta")),
0027       nintEta_(pset.getParameter<int>("nintEta")),
0028       useFabsEta_(pset.getParameter<bool>("useFabsEta")),
0029 
0030       //parameters for energy
0031       minEne_(pset.getParameter<double>("minEne")),
0032       maxEne_(pset.getParameter<double>("maxEne")),
0033       nintEne_(pset.getParameter<int>("nintEne")),
0034 
0035       //parameters for pt
0036       minPt_(pset.getParameter<double>("minPt")),
0037       maxPt_(pset.getParameter<double>("maxPt")),
0038       nintPt_(pset.getParameter<int>("nintPt")),
0039 
0040       //parameters for phi
0041       minPhi_(pset.getParameter<double>("minPhi")),
0042       maxPhi_(pset.getParameter<double>("maxPhi")),
0043       nintPhi_(pset.getParameter<int>("nintPhi")),
0044 
0045       //parameters for counting mixed hits SimClusters
0046       minMixedHitsSimCluster_(pset.getParameter<double>("minMixedHitsSimCluster")),
0047       maxMixedHitsSimCluster_(pset.getParameter<double>("maxMixedHitsSimCluster")),
0048       nintMixedHitsSimCluster_(pset.getParameter<int>("nintMixedHitsSimCluster")),
0049 
0050       //parameters for counting mixed hits clusters
0051       minMixedHitsCluster_(pset.getParameter<double>("minMixedHitsCluster")),
0052       maxMixedHitsCluster_(pset.getParameter<double>("maxMixedHitsCluster")),
0053       nintMixedHitsCluster_(pset.getParameter<int>("nintMixedHitsCluster")),
0054 
0055       //parameters for the total amount of energy clustered by all layer clusters (fraction over CaloParticless)
0056       minEneCl_(pset.getParameter<double>("minEneCl")),
0057       maxEneCl_(pset.getParameter<double>("maxEneCl")),
0058       nintEneCl_(pset.getParameter<int>("nintEneCl")),
0059 
0060       //parameters for the longitudinal depth barycenter.
0061       minLongDepBary_(pset.getParameter<double>("minLongDepBary")),
0062       maxLongDepBary_(pset.getParameter<double>("maxLongDepBary")),
0063       nintLongDepBary_(pset.getParameter<int>("nintLongDepBary")),
0064 
0065       //parameters for z positionof vertex plots
0066       minZpos_(pset.getParameter<double>("minZpos")),
0067       maxZpos_(pset.getParameter<double>("maxZpos")),
0068       nintZpos_(pset.getParameter<int>("nintZpos")),
0069 
0070       //Parameters for the total number of SimClusters per layer
0071       minTotNsimClsperlay_(pset.getParameter<double>("minTotNsimClsperlay")),
0072       maxTotNsimClsperlay_(pset.getParameter<double>("maxTotNsimClsperlay")),
0073       nintTotNsimClsperlay_(pset.getParameter<int>("nintTotNsimClsperlay")),
0074 
0075       //Parameters for the total number of layer clusters per layer
0076       minTotNClsperlay_(pset.getParameter<double>("minTotNClsperlay")),
0077       maxTotNClsperlay_(pset.getParameter<double>("maxTotNClsperlay")),
0078       nintTotNClsperlay_(pset.getParameter<int>("nintTotNClsperlay")),
0079 
0080       //Parameters for the energy clustered by layer clusters per layer (fraction over CaloParticless)
0081       minEneClperlay_(pset.getParameter<double>("minEneClperlay")),
0082       maxEneClperlay_(pset.getParameter<double>("maxEneClperlay")),
0083       nintEneClperlay_(pset.getParameter<int>("nintEneClperlay")),
0084 
0085       //Parameters for the score both for:
0086       //1. calo particle to layer clusters association per layer
0087       //2. layer cluster to calo particles association per layer
0088       minScore_(pset.getParameter<double>("minScore")),
0089       maxScore_(pset.getParameter<double>("maxScore")),
0090       nintScore_(pset.getParameter<int>("nintScore")),
0091 
0092       //Parameters for shared energy fraction. That is:
0093       //1. Fraction of each of the layer clusters energy related to a
0094       //calo particle over that calo particle's energy.
0095       //2. Fraction of each of the calo particles energy
0096       //related to a layer cluster over that layer cluster's energy.
0097       minSharedEneFrac_(pset.getParameter<double>("minSharedEneFrac")),
0098       maxSharedEneFrac_(pset.getParameter<double>("maxSharedEneFrac")),
0099       nintSharedEneFrac_(pset.getParameter<int>("nintSharedEneFrac")),
0100       minTSTSharedEneFracEfficiency_(pset.getParameter<double>("minTSTSharedEneFracEfficiency")),
0101 
0102       //Same as above for Tracksters
0103       minTSTSharedEneFrac_(pset.getParameter<double>("minTSTSharedEneFrac")),
0104       maxTSTSharedEneFrac_(pset.getParameter<double>("maxTSTSharedEneFrac")),
0105       nintTSTSharedEneFrac_(pset.getParameter<int>("nintTSTSharedEneFrac")),
0106 
0107       //Parameters for the total number of SimClusters per thickness
0108       minTotNsimClsperthick_(pset.getParameter<double>("minTotNsimClsperthick")),
0109       maxTotNsimClsperthick_(pset.getParameter<double>("maxTotNsimClsperthick")),
0110       nintTotNsimClsperthick_(pset.getParameter<int>("nintTotNsimClsperthick")),
0111 
0112       //Parameters for the total number of layer clusters per thickness
0113       minTotNClsperthick_(pset.getParameter<double>("minTotNClsperthick")),
0114       maxTotNClsperthick_(pset.getParameter<double>("maxTotNClsperthick")),
0115       nintTotNClsperthick_(pset.getParameter<int>("nintTotNClsperthick")),
0116 
0117       //Parameters for the total number of cells per per thickness per layer
0118       minTotNcellsperthickperlayer_(pset.getParameter<double>("minTotNcellsperthickperlayer")),
0119       maxTotNcellsperthickperlayer_(pset.getParameter<double>("maxTotNcellsperthickperlayer")),
0120       nintTotNcellsperthickperlayer_(pset.getParameter<int>("nintTotNcellsperthickperlayer")),
0121 
0122       //Parameters for the distance of cluster cells to seed cell per thickness per layer
0123       minDisToSeedperthickperlayer_(pset.getParameter<double>("minDisToSeedperthickperlayer")),
0124       maxDisToSeedperthickperlayer_(pset.getParameter<double>("maxDisToSeedperthickperlayer")),
0125       nintDisToSeedperthickperlayer_(pset.getParameter<int>("nintDisToSeedperthickperlayer")),
0126 
0127       //Parameters for the energy weighted distance of cluster cells to seed cell per thickness per layer
0128       minDisToSeedperthickperlayerenewei_(pset.getParameter<double>("minDisToSeedperthickperlayerenewei")),
0129       maxDisToSeedperthickperlayerenewei_(pset.getParameter<double>("maxDisToSeedperthickperlayerenewei")),
0130       nintDisToSeedperthickperlayerenewei_(pset.getParameter<int>("nintDisToSeedperthickperlayerenewei")),
0131 
0132       //Parameters for the distance of cluster cells to max cell per thickness per layer
0133       minDisToMaxperthickperlayer_(pset.getParameter<double>("minDisToMaxperthickperlayer")),
0134       maxDisToMaxperthickperlayer_(pset.getParameter<double>("maxDisToMaxperthickperlayer")),
0135       nintDisToMaxperthickperlayer_(pset.getParameter<int>("nintDisToMaxperthickperlayer")),
0136 
0137       //Parameters for the energy weighted distance of cluster cells to max cell per thickness per layer
0138       minDisToMaxperthickperlayerenewei_(pset.getParameter<double>("minDisToMaxperthickperlayerenewei")),
0139       maxDisToMaxperthickperlayerenewei_(pset.getParameter<double>("maxDisToMaxperthickperlayerenewei")),
0140       nintDisToMaxperthickperlayerenewei_(pset.getParameter<int>("nintDisToMaxperthickperlayerenewei")),
0141 
0142       //Parameters for the distance of seed cell to max cell per thickness per layer
0143       minDisSeedToMaxperthickperlayer_(pset.getParameter<double>("minDisSeedToMaxperthickperlayer")),
0144       maxDisSeedToMaxperthickperlayer_(pset.getParameter<double>("maxDisSeedToMaxperthickperlayer")),
0145       nintDisSeedToMaxperthickperlayer_(pset.getParameter<int>("nintDisSeedToMaxperthickperlayer")),
0146 
0147       //Parameters for the energy of a cluster per thickness per layer
0148       minClEneperthickperlayer_(pset.getParameter<double>("minClEneperthickperlayer")),
0149       maxClEneperthickperlayer_(pset.getParameter<double>("maxClEneperthickperlayer")),
0150       nintClEneperthickperlayer_(pset.getParameter<int>("nintClEneperthickperlayer")),
0151 
0152       //Parameters for the energy density of cluster cells per thickness
0153       minCellsEneDensperthick_(pset.getParameter<double>("minCellsEneDensperthick")),
0154       maxCellsEneDensperthick_(pset.getParameter<double>("maxCellsEneDensperthick")),
0155       nintCellsEneDensperthick_(pset.getParameter<int>("nintCellsEneDensperthick")),
0156 
0157       //Parameters for the total number of Tracksters per event
0158       // Always treat one event as two events, one in +z one in -z
0159       minTotNTSTs_(pset.getParameter<double>("minTotNTSTs")),
0160       maxTotNTSTs_(pset.getParameter<double>("maxTotNTSTs")),
0161       nintTotNTSTs_(pset.getParameter<int>("nintTotNTSTs")),
0162 
0163       //Parameters for the total number of layer clusters in Trackster
0164       minTotNClsinTSTs_(pset.getParameter<double>("minTotNClsinTSTs")),
0165       maxTotNClsinTSTs_(pset.getParameter<double>("maxTotNClsinTSTs")),
0166       nintTotNClsinTSTs_(pset.getParameter<int>("nintTotNClsinTSTs")),
0167 
0168       //Parameters for the total number of layer clusters in Trackster per layer
0169       minTotNClsinTSTsperlayer_(pset.getParameter<double>("minTotNClsinTSTsperlayer")),
0170       maxTotNClsinTSTsperlayer_(pset.getParameter<double>("maxTotNClsinTSTsperlayer")),
0171       nintTotNClsinTSTsperlayer_(pset.getParameter<int>("nintTotNClsinTSTsperlayer")),
0172 
0173       //Parameters for the multiplicity of layer clusters in Trackster
0174       minMplofLCs_(pset.getParameter<double>("minMplofLCs")),
0175       maxMplofLCs_(pset.getParameter<double>("maxMplofLCs")),
0176       nintMplofLCs_(pset.getParameter<int>("nintMplofLCs")),
0177 
0178       //Parameters for cluster size
0179       minSizeCLsinTSTs_(pset.getParameter<double>("minSizeCLsinTSTs")),
0180       maxSizeCLsinTSTs_(pset.getParameter<double>("maxSizeCLsinTSTs")),
0181       nintSizeCLsinTSTs_(pset.getParameter<int>("nintSizeCLsinTSTs")),
0182 
0183       //Parameters for the energy of a cluster per thickness per layer
0184       minClEnepermultiplicity_(pset.getParameter<double>("minClEnepermultiplicity")),
0185       maxClEnepermultiplicity_(pset.getParameter<double>("maxClEnepermultiplicity")),
0186       nintClEnepermultiplicity_(pset.getParameter<int>("nintClEnepermultiplicity")),
0187 
0188       //parameters for x
0189       minX_(pset.getParameter<double>("minX")),
0190       maxX_(pset.getParameter<double>("maxX")),
0191       nintX_(pset.getParameter<int>("nintX")),
0192 
0193       //parameters for y
0194       minY_(pset.getParameter<double>("minY")),
0195       maxY_(pset.getParameter<double>("maxY")),
0196       nintY_(pset.getParameter<int>("nintY")),
0197 
0198       //parameters for z
0199       minZ_(pset.getParameter<double>("minZ")),
0200       maxZ_(pset.getParameter<double>("maxZ")),
0201       nintZ_(pset.getParameter<int>("nintZ")) {}
0202 
0203 HGVHistoProducerAlgo::~HGVHistoProducerAlgo() {}
0204 
0205 void HGVHistoProducerAlgo::bookInfo(DQMStore::IBooker& ibook, Histograms& histograms) {
0206   histograms.lastLayerEEzm = ibook.bookInt("lastLayerEEzm");
0207   histograms.lastLayerFHzm = ibook.bookInt("lastLayerFHzm");
0208   histograms.maxlayerzm = ibook.bookInt("maxlayerzm");
0209   histograms.lastLayerEEzp = ibook.bookInt("lastLayerEEzp");
0210   histograms.lastLayerFHzp = ibook.bookInt("lastLayerFHzp");
0211   histograms.maxlayerzp = ibook.bookInt("maxlayerzp");
0212 }
0213 
0214 void HGVHistoProducerAlgo::bookCaloParticleHistos(DQMStore::IBooker& ibook,
0215                                                   Histograms& histograms,
0216                                                   int pdgid,
0217                                                   unsigned int layers) {
0218   histograms.h_caloparticle_eta[pdgid] =
0219       ibook.book1D("N of caloparticle vs eta", "N of caloParticles vs eta", nintEta_, minEta_, maxEta_);
0220   histograms.h_caloparticle_eta_Zorigin[pdgid] =
0221       ibook.book2D("Eta vs Zorigin", "Eta vs Zorigin", nintEta_, minEta_, maxEta_, nintZpos_, minZpos_, maxZpos_);
0222 
0223   histograms.h_caloparticle_energy[pdgid] =
0224       ibook.book1D("Energy", "Energy of CaloParticles; Energy [GeV]", nintEne_, minEne_, maxEne_);
0225   histograms.h_caloparticle_pt[pdgid] = ibook.book1D("Pt", "Pt of CaloParticles", nintPt_, minPt_, maxPt_);
0226   histograms.h_caloparticle_phi[pdgid] = ibook.book1D("Phi", "Phi of CaloParticles", nintPhi_, minPhi_, maxPhi_);
0227   histograms.h_caloparticle_selfenergy[pdgid] =
0228       ibook.book1D("SelfEnergy", "Total Energy of Hits in Sim Clusters (matched)", nintEne_, minEne_, maxEne_);
0229   histograms.h_caloparticle_energyDifference[pdgid] =
0230       ibook.book1D("EnergyDifference", "(Energy-SelfEnergy)/Energy", 300., -5., 1.);
0231 
0232   histograms.h_caloparticle_nSimClusters[pdgid] =
0233       ibook.book1D("Num Sim Clusters", "Num Sim Clusters in CaloParticles", 100, 0., 100.);
0234   histograms.h_caloparticle_nHitsInSimClusters[pdgid] =
0235       ibook.book1D("Num Hits in Sim Clusters", "Num Hits in Sim Clusters in CaloParticles", 1000, 0., 1000.);
0236   histograms.h_caloparticle_nHitsInSimClusters_matchedtoRecHit[pdgid] = ibook.book1D(
0237       "Num Rec-matched Hits in Sim Clusters", "Num Hits in Sim Clusters (matched) in CaloParticles", 1000, 0., 1000.);
0238 
0239   histograms.h_caloparticle_nHits_matched_energy[pdgid] =
0240       ibook.book1D("Energy of Rec-matched Hits", "Energy of Hits in Sim Clusters (matched)", 100, 0., 10.);
0241   histograms.h_caloparticle_nHits_matched_energy_layer[pdgid] =
0242       ibook.book2D("Energy of Rec-matched Hits vs layer",
0243                    "Energy of Hits in Sim Clusters (matched) vs layer",
0244                    2 * layers,
0245                    0.,
0246                    (float)2 * layers,
0247                    100,
0248                    0.,
0249                    10.);
0250   histograms.h_caloparticle_nHits_matched_energy_layer_1SimCl[pdgid] =
0251       ibook.book2D("Energy of Rec-matched Hits vs layer (1SC)",
0252                    "Energy of Hits only 1 Sim Clusters (matched) vs layer",
0253                    2 * layers,
0254                    0.,
0255                    (float)2 * layers,
0256                    100,
0257                    0.,
0258                    10.);
0259   histograms.h_caloparticle_sum_energy_layer[pdgid] =
0260       ibook.book2D("Rec-matched Hits Sum Energy vs layer",
0261                    "Rescaled Sum Energy of Hits in Sim Clusters (matched) vs layer",
0262                    2 * layers,
0263                    0.,
0264                    (float)2 * layers,
0265                    110,
0266                    0.,
0267                    110.);
0268   histograms.h_caloparticle_fractions[pdgid] =
0269       ibook.book2D("HitFractions", "Hit fractions;Hit fraction;E_{hit}^{2} fraction", 101, 0, 1.01, 100, 0, 1);
0270   histograms.h_caloparticle_fractions_weight[pdgid] = ibook.book2D(
0271       "HitFractions_weighted", "Hit fractions weighted;Hit fraction;E_{hit}^{2} fraction", 101, 0, 1.01, 100, 0, 1);
0272 
0273   histograms.h_caloparticle_firstlayer[pdgid] =
0274       ibook.book1D("First Layer", "First layer of the CaloParticles", 2 * layers, 0., (float)2 * layers);
0275   histograms.h_caloparticle_lastlayer[pdgid] =
0276       ibook.book1D("Last Layer", "Last layer of the CaloParticles", 2 * layers, 0., (float)2 * layers);
0277   histograms.h_caloparticle_layersnum[pdgid] =
0278       ibook.book1D("Number of Layers", "Number of layers of the CaloParticles", 2 * layers, 0., (float)2 * layers);
0279   histograms.h_caloparticle_firstlayer_matchedtoRecHit[pdgid] = ibook.book1D(
0280       "First Layer (rec-matched hit)", "First layer of the CaloParticles (matched)", 2 * layers, 0., (float)2 * layers);
0281   histograms.h_caloparticle_lastlayer_matchedtoRecHit[pdgid] = ibook.book1D(
0282       "Last Layer (rec-matched hit)", "Last layer of the CaloParticles (matched)", 2 * layers, 0., (float)2 * layers);
0283   histograms.h_caloparticle_layersnum_matchedtoRecHit[pdgid] =
0284       ibook.book1D("Number of Layers (rec-matched hit)",
0285                    "Number of layers of the CaloParticles (matched)",
0286                    2 * layers,
0287                    0.,
0288                    (float)2 * layers);
0289 }
0290 
0291 void HGVHistoProducerAlgo::bookSimClusterHistos(DQMStore::IBooker& ibook,
0292                                                 Histograms& histograms,
0293                                                 unsigned int layers,
0294                                                 std::vector<int> thicknesses) {
0295   //---------------------------------------------------------------------------------------------------------------------------
0296   for (unsigned ilayer = 0; ilayer < 2 * layers; ++ilayer) {
0297     auto istr1 = std::to_string(ilayer);
0298     while (istr1.size() < 2) {
0299       istr1.insert(0, "0");
0300     }
0301     // Make a mapping to the regural layer naming plus z- or z+ for convenience
0302     std::string istr2 = "";
0303     // first with the -z endcap
0304     if (ilayer < layers) {
0305       istr2 = std::to_string(ilayer + 1) + " in z-";
0306     } else {  // then for the +z
0307       istr2 = std::to_string(ilayer - (layers - 1)) + " in z+";
0308     }
0309     histograms.h_simclusternum_perlayer[ilayer] = ibook.book1D("totsimclusternum_layer_" + istr1,
0310                                                                "total number of SimClusters for layer " + istr2,
0311                                                                nintTotNsimClsperlay_,
0312                                                                minTotNsimClsperlay_,
0313                                                                maxTotNsimClsperlay_);
0314 
0315   }  //end of loop over layers
0316   //---------------------------------------------------------------------------------------------------------------------------
0317   for (std::vector<int>::iterator it = thicknesses.begin(); it != thicknesses.end(); ++it) {
0318     auto istr = std::to_string(*it);
0319     histograms.h_simclusternum_perthick[(*it)] = ibook.book1D("totsimclusternum_thick_" + istr,
0320                                                               "total number of simclusters for thickness " + istr,
0321                                                               nintTotNsimClsperthick_,
0322                                                               minTotNsimClsperthick_,
0323                                                               maxTotNsimClsperthick_);
0324   }  //end of loop over thicknesses
0325 
0326   //---------------------------------------------------------------------------------------------------------------------------
0327   //z-
0328   histograms.h_mixedhitssimcluster_zminus =
0329       ibook.book1D("mixedhitssimcluster_zminus",
0330                    "N of simclusters that contain hits of more than one kind in z-",
0331                    nintMixedHitsSimCluster_,
0332                    minMixedHitsSimCluster_,
0333                    maxMixedHitsSimCluster_);
0334   //z+
0335   histograms.h_mixedhitssimcluster_zplus =
0336       ibook.book1D("mixedhitssimcluster_zplus",
0337                    "N of simclusters that contain hits of more than one kind in z+",
0338                    nintMixedHitsSimCluster_,
0339                    minMixedHitsSimCluster_,
0340                    maxMixedHitsSimCluster_);
0341 }
0342 
0343 void HGVHistoProducerAlgo::bookSimClusterAssociationHistos(DQMStore::IBooker& ibook,
0344                                                            Histograms& histograms,
0345                                                            unsigned int layers,
0346                                                            std::vector<int> thicknesses) {
0347   std::unordered_map<int, dqm::reco::MonitorElement*> denom_layercl_in_simcl_eta_perlayer;
0348   denom_layercl_in_simcl_eta_perlayer.clear();
0349   std::unordered_map<int, dqm::reco::MonitorElement*> denom_layercl_in_simcl_phi_perlayer;
0350   denom_layercl_in_simcl_phi_perlayer.clear();
0351   std::unordered_map<int, dqm::reco::MonitorElement*> score_layercl2simcluster_perlayer;
0352   score_layercl2simcluster_perlayer.clear();
0353   std::unordered_map<int, dqm::reco::MonitorElement*> sharedenergy_layercl2simcluster_perlayer;
0354   sharedenergy_layercl2simcluster_perlayer.clear();
0355   std::unordered_map<int, dqm::reco::MonitorElement*> energy_vs_score_layercl2simcluster_perlayer;
0356   energy_vs_score_layercl2simcluster_perlayer.clear();
0357   std::unordered_map<int, dqm::reco::MonitorElement*> num_layercl_in_simcl_eta_perlayer;
0358   num_layercl_in_simcl_eta_perlayer.clear();
0359   std::unordered_map<int, dqm::reco::MonitorElement*> num_layercl_in_simcl_phi_perlayer;
0360   num_layercl_in_simcl_phi_perlayer.clear();
0361   std::unordered_map<int, dqm::reco::MonitorElement*> numMerge_layercl_in_simcl_eta_perlayer;
0362   numMerge_layercl_in_simcl_eta_perlayer.clear();
0363   std::unordered_map<int, dqm::reco::MonitorElement*> numMerge_layercl_in_simcl_phi_perlayer;
0364   numMerge_layercl_in_simcl_phi_perlayer.clear();
0365   std::unordered_map<int, dqm::reco::MonitorElement*> sharedenergy_layercl2simcluster_vs_eta_perlayer;
0366   sharedenergy_layercl2simcluster_vs_eta_perlayer.clear();
0367   std::unordered_map<int, dqm::reco::MonitorElement*> sharedenergy_layercl2simcluster_vs_phi_perlayer;
0368   sharedenergy_layercl2simcluster_vs_phi_perlayer.clear();
0369   std::unordered_map<int, dqm::reco::MonitorElement*> denom_simcluster_eta_perlayer;
0370   denom_simcluster_eta_perlayer.clear();
0371   std::unordered_map<int, dqm::reco::MonitorElement*> denom_simcluster_phi_perlayer;
0372   denom_simcluster_phi_perlayer.clear();
0373   std::unordered_map<int, dqm::reco::MonitorElement*> score_simcluster2layercl_perlayer;
0374   score_simcluster2layercl_perlayer.clear();
0375   std::unordered_map<int, dqm::reco::MonitorElement*> sharedenergy_simcluster2layercl_perlayer;
0376   sharedenergy_simcluster2layercl_perlayer.clear();
0377   std::unordered_map<int, dqm::reco::MonitorElement*> energy_vs_score_simcluster2layercl_perlayer;
0378   energy_vs_score_simcluster2layercl_perlayer.clear();
0379   std::unordered_map<int, dqm::reco::MonitorElement*> num_simcluster_eta_perlayer;
0380   num_simcluster_eta_perlayer.clear();
0381   std::unordered_map<int, dqm::reco::MonitorElement*> num_simcluster_phi_perlayer;
0382   num_simcluster_phi_perlayer.clear();
0383   std::unordered_map<int, dqm::reco::MonitorElement*> numDup_simcluster_eta_perlayer;
0384   numDup_simcluster_eta_perlayer.clear();
0385   std::unordered_map<int, dqm::reco::MonitorElement*> numDup_simcluster_phi_perlayer;
0386   numDup_simcluster_phi_perlayer.clear();
0387   std::unordered_map<int, dqm::reco::MonitorElement*> sharedenergy_simcluster2layercl_vs_eta_perlayer;
0388   sharedenergy_simcluster2layercl_vs_eta_perlayer.clear();
0389   std::unordered_map<int, dqm::reco::MonitorElement*> sharedenergy_simcluster2layercl_vs_phi_perlayer;
0390   sharedenergy_simcluster2layercl_vs_phi_perlayer.clear();
0391 
0392   //---------------------------------------------------------------------------------------------------------------------------
0393   for (unsigned ilayer = 0; ilayer < 2 * layers; ++ilayer) {
0394     auto istr1 = std::to_string(ilayer);
0395     while (istr1.size() < 2) {
0396       istr1.insert(0, "0");
0397     }
0398     // Make a mapping to the regural layer naming plus z- or z+ for convenience
0399     std::string istr2 = "";
0400     // first with the -z endcap
0401     if (ilayer < layers) {
0402       istr2 = std::to_string(ilayer + 1) + " in z-";
0403     } else {  // then for the +z
0404       istr2 = std::to_string(ilayer - (layers - 1)) + " in z+";
0405     }
0406     //-------------------------------------------------------------------------------------------------------------------------
0407     denom_layercl_in_simcl_eta_perlayer[ilayer] =
0408         ibook.book1D("Denom_LayerCluster_in_SimCluster_Eta_perlayer" + istr1,
0409                      "Denom LayerCluster in SimCluster Eta per Layer Cluster for layer " + istr2,
0410                      nintEta_,
0411                      minEta_,
0412                      maxEta_);
0413     //-------------------------------------------------------------------------------------------------------------------------
0414     denom_layercl_in_simcl_phi_perlayer[ilayer] =
0415         ibook.book1D("Denom_LayerCluster_in_SimCluster_Phi_perlayer" + istr1,
0416                      "Denom LayerCluster in SimCluster Phi per Layer Cluster for layer " + istr2,
0417                      nintPhi_,
0418                      minPhi_,
0419                      maxPhi_);
0420     //-------------------------------------------------------------------------------------------------------------------------
0421     score_layercl2simcluster_perlayer[ilayer] = ibook.book1D("Score_layercl2simcluster_perlayer" + istr1,
0422                                                              "Score of Layer Cluster per SimCluster for layer " + istr2,
0423                                                              nintScore_,
0424                                                              minScore_,
0425                                                              maxScore_);
0426     //-------------------------------------------------------------------------------------------------------------------------
0427     score_simcluster2layercl_perlayer[ilayer] = ibook.book1D("Score_simcluster2layercl_perlayer" + istr1,
0428                                                              "Score of SimCluster per Layer Cluster for layer " + istr2,
0429                                                              nintScore_,
0430                                                              minScore_,
0431                                                              maxScore_);
0432     //-------------------------------------------------------------------------------------------------------------------------
0433     energy_vs_score_simcluster2layercl_perlayer[ilayer] =
0434         ibook.book2D("Energy_vs_Score_simcluster2layer_perlayer" + istr1,
0435                      "Energy vs Score of SimCluster per Layer Cluster for layer " + istr2,
0436                      nintScore_,
0437                      minScore_,
0438                      maxScore_,
0439                      nintSharedEneFrac_,
0440                      minSharedEneFrac_,
0441                      maxSharedEneFrac_);
0442     //-------------------------------------------------------------------------------------------------------------------------
0443     energy_vs_score_layercl2simcluster_perlayer[ilayer] =
0444         ibook.book2D("Energy_vs_Score_layer2simcluster_perlayer" + istr1,
0445                      "Energy vs Score of Layer Cluster per SimCluster for layer " + istr2,
0446                      nintScore_,
0447                      minScore_,
0448                      maxScore_,
0449                      nintSharedEneFrac_,
0450                      minSharedEneFrac_,
0451                      maxSharedEneFrac_);
0452     //-------------------------------------------------------------------------------------------------------------------------
0453     sharedenergy_simcluster2layercl_perlayer[ilayer] =
0454         ibook.book1D("SharedEnergy_simcluster2layercl_perlayer" + istr1,
0455                      "Shared Energy of SimCluster per Layer Cluster for layer " + istr2,
0456                      nintSharedEneFrac_,
0457                      minSharedEneFrac_,
0458                      maxSharedEneFrac_);
0459     //-------------------------------------------------------------------------------------------------------------------------
0460     sharedenergy_simcluster2layercl_vs_eta_perlayer[ilayer] =
0461         ibook.bookProfile("SharedEnergy_simcluster2layercl_vs_eta_perlayer" + istr1,
0462                           "Shared Energy of SimCluster vs #eta per best Layer Cluster for layer " + istr2,
0463                           nintEta_,
0464                           minEta_,
0465                           maxEta_,
0466                           minSharedEneFrac_,
0467                           maxSharedEneFrac_);
0468     //-------------------------------------------------------------------------------------------------------------------------
0469     sharedenergy_simcluster2layercl_vs_phi_perlayer[ilayer] =
0470         ibook.bookProfile("SharedEnergy_simcluster2layercl_vs_phi_perlayer" + istr1,
0471                           "Shared Energy of SimCluster vs #phi per best Layer Cluster for layer " + istr2,
0472                           nintPhi_,
0473                           minPhi_,
0474                           maxPhi_,
0475                           minSharedEneFrac_,
0476                           maxSharedEneFrac_);
0477     //-------------------------------------------------------------------------------------------------------------------------
0478     sharedenergy_layercl2simcluster_perlayer[ilayer] =
0479         ibook.book1D("SharedEnergy_layercluster2simcluster_perlayer" + istr1,
0480                      "Shared Energy of Layer Cluster per SimCluster for layer " + istr2,
0481                      nintSharedEneFrac_,
0482                      minSharedEneFrac_,
0483                      maxSharedEneFrac_);
0484     //-------------------------------------------------------------------------------------------------------------------------
0485     sharedenergy_layercl2simcluster_vs_eta_perlayer[ilayer] =
0486         ibook.bookProfile("SharedEnergy_layercl2simcluster_vs_eta_perlayer" + istr1,
0487                           "Shared Energy of LayerCluster vs #eta per best SimCluster for layer " + istr2,
0488                           nintEta_,
0489                           minEta_,
0490                           maxEta_,
0491                           minSharedEneFrac_,
0492                           maxSharedEneFrac_);
0493     //-------------------------------------------------------------------------------------------------------------------------
0494     sharedenergy_layercl2simcluster_vs_phi_perlayer[ilayer] =
0495         ibook.bookProfile("SharedEnergy_layercl2simcluster_vs_phi_perlayer" + istr1,
0496                           "Shared Energy of LayerCluster vs #phi per best SimCluster for layer " + istr2,
0497                           nintPhi_,
0498                           minPhi_,
0499                           maxPhi_,
0500                           minSharedEneFrac_,
0501                           maxSharedEneFrac_);
0502     //-------------------------------------------------------------------------------------------------------------------------
0503     num_simcluster_eta_perlayer[ilayer] = ibook.book1D("Num_SimCluster_Eta_perlayer" + istr1,
0504                                                        "Num SimCluster Eta per Layer Cluster for layer " + istr2,
0505                                                        nintEta_,
0506                                                        minEta_,
0507                                                        maxEta_);
0508     //-------------------------------------------------------------------------------------------------------------------------
0509     numDup_simcluster_eta_perlayer[ilayer] =
0510         ibook.book1D("NumDup_SimCluster_Eta_perlayer" + istr1,
0511                      "Num Duplicate SimCluster Eta per Layer Cluster for layer " + istr2,
0512                      nintEta_,
0513                      minEta_,
0514                      maxEta_);
0515     //-------------------------------------------------------------------------------------------------------------------------
0516     denom_simcluster_eta_perlayer[ilayer] = ibook.book1D("Denom_SimCluster_Eta_perlayer" + istr1,
0517                                                          "Denom SimCluster Eta per Layer Cluster for layer " + istr2,
0518                                                          nintEta_,
0519                                                          minEta_,
0520                                                          maxEta_);
0521     //-------------------------------------------------------------------------------------------------------------------------
0522     num_simcluster_phi_perlayer[ilayer] = ibook.book1D("Num_SimCluster_Phi_perlayer" + istr1,
0523                                                        "Num SimCluster Phi per Layer Cluster for layer " + istr2,
0524                                                        nintPhi_,
0525                                                        minPhi_,
0526                                                        maxPhi_);
0527     //-------------------------------------------------------------------------------------------------------------------------
0528     numDup_simcluster_phi_perlayer[ilayer] =
0529         ibook.book1D("NumDup_SimCluster_Phi_perlayer" + istr1,
0530                      "Num Duplicate SimCluster Phi per Layer Cluster for layer " + istr2,
0531                      nintPhi_,
0532                      minPhi_,
0533                      maxPhi_);
0534     //-------------------------------------------------------------------------------------------------------------------------
0535     denom_simcluster_phi_perlayer[ilayer] = ibook.book1D("Denom_SimCluster_Phi_perlayer" + istr1,
0536                                                          "Denom SimCluster Phi per Layer Cluster for layer " + istr2,
0537                                                          nintPhi_,
0538                                                          minPhi_,
0539                                                          maxPhi_);
0540     //-------------------------------------------------------------------------------------------------------------------------
0541     num_layercl_in_simcl_eta_perlayer[ilayer] =
0542         ibook.book1D("Num_LayerCluster_in_SimCluster_Eta_perlayer" + istr1,
0543                      "Num LayerCluster Eta per Layer Cluster in SimCluster for layer " + istr2,
0544                      nintEta_,
0545                      minEta_,
0546                      maxEta_);
0547     //-------------------------------------------------------------------------------------------------------------------------
0548     numMerge_layercl_in_simcl_eta_perlayer[ilayer] =
0549         ibook.book1D("NumMerge_LayerCluster_in_SimCluster_Eta_perlayer" + istr1,
0550                      "Num Merge LayerCluster Eta per Layer Cluster in SimCluster for layer " + istr2,
0551                      nintEta_,
0552                      minEta_,
0553                      maxEta_);
0554     //-------------------------------------------------------------------------------------------------------------------------
0555     num_layercl_in_simcl_phi_perlayer[ilayer] =
0556         ibook.book1D("Num_LayerCluster_in_SimCluster_Phi_perlayer" + istr1,
0557                      "Num LayerCluster Phi per Layer Cluster in SimCluster for layer " + istr2,
0558                      nintPhi_,
0559                      minPhi_,
0560                      maxPhi_);
0561     //-------------------------------------------------------------------------------------------------------------------------
0562     numMerge_layercl_in_simcl_phi_perlayer[ilayer] =
0563         ibook.book1D("NumMerge_LayerCluster_in_SimCluster_Phi_perlayer" + istr1,
0564                      "Num Merge LayerCluster Phi per Layer Cluster in SimCluster for layer " + istr2,
0565                      nintPhi_,
0566                      minPhi_,
0567                      maxPhi_);
0568 
0569   }  //end of loop over layers
0570 
0571   histograms.h_denom_layercl_in_simcl_eta_perlayer.push_back(std::move(denom_layercl_in_simcl_eta_perlayer));
0572   histograms.h_denom_layercl_in_simcl_phi_perlayer.push_back(std::move(denom_layercl_in_simcl_phi_perlayer));
0573   histograms.h_score_layercl2simcluster_perlayer.push_back(std::move(score_layercl2simcluster_perlayer));
0574   histograms.h_sharedenergy_layercl2simcluster_perlayer.push_back(std::move(sharedenergy_layercl2simcluster_perlayer));
0575   histograms.h_energy_vs_score_layercl2simcluster_perlayer.push_back(
0576       std::move(energy_vs_score_layercl2simcluster_perlayer));
0577   histograms.h_num_layercl_in_simcl_eta_perlayer.push_back(std::move(num_layercl_in_simcl_eta_perlayer));
0578   histograms.h_num_layercl_in_simcl_phi_perlayer.push_back(std::move(num_layercl_in_simcl_phi_perlayer));
0579   histograms.h_numMerge_layercl_in_simcl_eta_perlayer.push_back(std::move(numMerge_layercl_in_simcl_eta_perlayer));
0580   histograms.h_numMerge_layercl_in_simcl_phi_perlayer.push_back(std::move(numMerge_layercl_in_simcl_phi_perlayer));
0581   histograms.h_sharedenergy_layercl2simcluster_vs_eta_perlayer.push_back(
0582       std::move(sharedenergy_layercl2simcluster_vs_eta_perlayer));
0583   histograms.h_sharedenergy_layercl2simcluster_vs_phi_perlayer.push_back(
0584       std::move(sharedenergy_layercl2simcluster_vs_phi_perlayer));
0585   histograms.h_denom_simcluster_eta_perlayer.push_back(std::move(denom_simcluster_eta_perlayer));
0586   histograms.h_denom_simcluster_phi_perlayer.push_back(std::move(denom_simcluster_phi_perlayer));
0587   histograms.h_score_simcluster2layercl_perlayer.push_back(std::move(score_simcluster2layercl_perlayer));
0588   histograms.h_sharedenergy_simcluster2layercl_perlayer.push_back(std::move(sharedenergy_simcluster2layercl_perlayer));
0589   histograms.h_energy_vs_score_simcluster2layercl_perlayer.push_back(
0590       std::move(energy_vs_score_simcluster2layercl_perlayer));
0591   histograms.h_num_simcluster_eta_perlayer.push_back(std::move(num_simcluster_eta_perlayer));
0592   histograms.h_num_simcluster_phi_perlayer.push_back(std::move(num_simcluster_phi_perlayer));
0593   histograms.h_numDup_simcluster_eta_perlayer.push_back(std::move(numDup_simcluster_eta_perlayer));
0594   histograms.h_numDup_simcluster_phi_perlayer.push_back(std::move(numDup_simcluster_phi_perlayer));
0595   histograms.h_sharedenergy_simcluster2layercl_vs_eta_perlayer.push_back(
0596       std::move(sharedenergy_simcluster2layercl_vs_eta_perlayer));
0597   histograms.h_sharedenergy_simcluster2layercl_vs_phi_perlayer.push_back(
0598       std::move(sharedenergy_simcluster2layercl_vs_phi_perlayer));
0599 }
0600 void HGVHistoProducerAlgo::bookClusterHistos_ClusterLevel(DQMStore::IBooker& ibook,
0601                                                           Histograms& histograms,
0602                                                           unsigned int layers,
0603                                                           std::vector<int> thicknesses,
0604                                                           std::string pathtomatbudfile) {
0605   //---------------------------------------------------------------------------------------------------------------------------
0606   histograms.h_cluster_eta.push_back(
0607       ibook.book1D("num_reco_cluster_eta", "N of reco clusters vs eta", nintEta_, minEta_, maxEta_));
0608 
0609   //---------------------------------------------------------------------------------------------------------------------------
0610   //z-
0611   histograms.h_mixedhitscluster_zminus.push_back(
0612       ibook.book1D("mixedhitscluster_zminus",
0613                    "N of reco clusters that contain hits of more than one kind in z-",
0614                    nintMixedHitsCluster_,
0615                    minMixedHitsCluster_,
0616                    maxMixedHitsCluster_));
0617   //z+
0618   histograms.h_mixedhitscluster_zplus.push_back(
0619       ibook.book1D("mixedhitscluster_zplus",
0620                    "N of reco clusters that contain hits of more than one kind in z+",
0621                    nintMixedHitsCluster_,
0622                    minMixedHitsCluster_,
0623                    maxMixedHitsCluster_));
0624 
0625   //---------------------------------------------------------------------------------------------------------------------------
0626   //z-
0627   histograms.h_energyclustered_zminus.push_back(
0628       ibook.book1D("energyclustered_zminus",
0629                    "percent of total energy clustered by all layer clusters over CaloParticless energy in z-",
0630                    nintEneCl_,
0631                    minEneCl_,
0632                    maxEneCl_));
0633   //z+
0634   histograms.h_energyclustered_zplus.push_back(
0635       ibook.book1D("energyclustered_zplus",
0636                    "percent of total energy clustered by all layer clusters over CaloParticless energy in z+",
0637                    nintEneCl_,
0638                    minEneCl_,
0639                    maxEneCl_));
0640 
0641   //---------------------------------------------------------------------------------------------------------------------------
0642   //z-
0643   std::string subpathtomat = pathtomatbudfile.substr(pathtomatbudfile.find("Validation"));
0644   histograms.h_longdepthbarycentre_zminus.push_back(
0645       ibook.book1D("longdepthbarycentre_zminus",
0646                    "The longitudinal depth barycentre in z- for " + subpathtomat,
0647                    nintLongDepBary_,
0648                    minLongDepBary_,
0649                    maxLongDepBary_));
0650   //z+
0651   histograms.h_longdepthbarycentre_zplus.push_back(
0652       ibook.book1D("longdepthbarycentre_zplus",
0653                    "The longitudinal depth barycentre in z+ for " + subpathtomat,
0654                    nintLongDepBary_,
0655                    minLongDepBary_,
0656                    maxLongDepBary_));
0657 
0658   //---------------------------------------------------------------------------------------------------------------------------
0659   for (unsigned ilayer = 0; ilayer < 2 * layers; ++ilayer) {
0660     auto istr1 = std::to_string(ilayer);
0661     while (istr1.size() < 2) {
0662       istr1.insert(0, "0");
0663     }
0664     // Make a mapping to the regural layer naming plus z- or z+ for convenience
0665     std::string istr2 = "";
0666     // first with the -z endcap
0667     if (ilayer < layers) {
0668       istr2 = std::to_string(ilayer + 1) + " in z-";
0669     } else {  // then for the +z
0670       istr2 = std::to_string(ilayer - (layers - 1)) + " in z+";
0671     }
0672     histograms.h_clusternum_perlayer[ilayer] = ibook.book1D("totclusternum_layer_" + istr1,
0673                                                             "total number of layer clusters for layer " + istr2,
0674                                                             nintTotNClsperlay_,
0675                                                             minTotNClsperlay_,
0676                                                             maxTotNClsperlay_);
0677     histograms.h_energyclustered_perlayer[ilayer] = ibook.book1D(
0678         "energyclustered_perlayer" + istr1,
0679         "percent of total energy clustered by layer clusters over CaloParticless energy for layer " + istr2,
0680         nintEneClperlay_,
0681         minEneClperlay_,
0682         maxEneClperlay_);
0683   }
0684 
0685   //---------------------------------------------------------------------------------------------------------------------------
0686   for (std::vector<int>::iterator it = thicknesses.begin(); it != thicknesses.end(); ++it) {
0687     auto istr = std::to_string(*it);
0688     histograms.h_clusternum_perthick[(*it)] = ibook.book1D("totclusternum_thick_" + istr,
0689                                                            "total number of layer clusters for thickness " + istr,
0690                                                            nintTotNClsperthick_,
0691                                                            minTotNClsperthick_,
0692                                                            maxTotNClsperthick_);
0693   }
0694   //---------------------------------------------------------------------------------------------------------------------------
0695 }
0696 
0697 void HGVHistoProducerAlgo::bookClusterHistos_LCtoCP_association(DQMStore::IBooker& ibook,
0698                                                                 Histograms& histograms,
0699                                                                 unsigned int layers) {
0700   //----------------------------------------------------------------------------------------------------------------------------
0701   for (unsigned ilayer = 0; ilayer < 2 * layers; ++ilayer) {
0702     auto istr1 = std::to_string(ilayer);
0703     while (istr1.size() < 2) {
0704       istr1.insert(0, "0");
0705     }
0706     // Make a mapping to the regural layer naming plus z- or z+ for convenience
0707     std::string istr2 = "";
0708     // first with the -z endcap
0709     if (ilayer < layers) {
0710       istr2 = std::to_string(ilayer + 1) + " in z-";
0711     } else {  // then for the +z
0712       istr2 = std::to_string(ilayer - (layers - 1)) + " in z+";
0713     }
0714     histograms.h_score_layercl2caloparticle_perlayer[ilayer] =
0715         ibook.book1D("Score_layercl2caloparticle_perlayer" + istr1,
0716                      "Score of Layer Cluster per CaloParticle for layer " + istr2,
0717                      nintScore_,
0718                      minScore_,
0719                      maxScore_);
0720     histograms.h_score_caloparticle2layercl_perlayer[ilayer] =
0721         ibook.book1D("Score_caloparticle2layercl_perlayer" + istr1,
0722                      "Score of CaloParticle per Layer Cluster for layer " + istr2,
0723                      nintScore_,
0724                      minScore_,
0725                      maxScore_);
0726     histograms.h_energy_vs_score_caloparticle2layercl_perlayer[ilayer] =
0727         ibook.book2D("Energy_vs_Score_caloparticle2layer_perlayer" + istr1,
0728                      "Energy vs Score of CaloParticle per Layer Cluster for layer " + istr2,
0729                      nintScore_,
0730                      minScore_,
0731                      maxScore_,
0732                      nintSharedEneFrac_,
0733                      minSharedEneFrac_,
0734                      maxSharedEneFrac_);
0735     histograms.h_energy_vs_score_layercl2caloparticle_perlayer[ilayer] =
0736         ibook.book2D("Energy_vs_Score_layer2caloparticle_perlayer" + istr1,
0737                      "Energy vs Score of Layer Cluster per CaloParticle Layer for layer " + istr2,
0738                      nintScore_,
0739                      minScore_,
0740                      maxScore_,
0741                      nintSharedEneFrac_,
0742                      minSharedEneFrac_,
0743                      maxSharedEneFrac_);
0744     histograms.h_sharedenergy_caloparticle2layercl_perlayer[ilayer] =
0745         ibook.book1D("SharedEnergy_caloparticle2layercl_perlayer" + istr1,
0746                      "Shared Energy of CaloParticle per Layer Cluster for layer " + istr2,
0747                      nintSharedEneFrac_,
0748                      minSharedEneFrac_,
0749                      maxSharedEneFrac_);
0750     histograms.h_sharedenergy_caloparticle2layercl_vs_eta_perlayer[ilayer] =
0751         ibook.bookProfile("SharedEnergy_caloparticle2layercl_vs_eta_perlayer" + istr1,
0752                           "Shared Energy of CaloParticle vs #eta per best Layer Cluster for layer " + istr2,
0753                           nintEta_,
0754                           minEta_,
0755                           maxEta_,
0756                           minSharedEneFrac_,
0757                           maxSharedEneFrac_);
0758     histograms.h_sharedenergy_caloparticle2layercl_vs_phi_perlayer[ilayer] =
0759         ibook.bookProfile("SharedEnergy_caloparticle2layercl_vs_phi_perlayer" + istr1,
0760                           "Shared Energy of CaloParticle vs #phi per best Layer Cluster for layer " + istr2,
0761                           nintPhi_,
0762                           minPhi_,
0763                           maxPhi_,
0764                           minSharedEneFrac_,
0765                           maxSharedEneFrac_);
0766     histograms.h_sharedenergy_layercl2caloparticle_perlayer[ilayer] =
0767         ibook.book1D("SharedEnergy_layercluster2caloparticle_perlayer" + istr1,
0768                      "Shared Energy of Layer Cluster per Layer Calo Particle for layer " + istr2,
0769                      nintSharedEneFrac_,
0770                      minSharedEneFrac_,
0771                      maxSharedEneFrac_);
0772     histograms.h_sharedenergy_layercl2caloparticle_vs_eta_perlayer[ilayer] =
0773         ibook.bookProfile("SharedEnergy_layercl2caloparticle_vs_eta_perlayer" + istr1,
0774                           "Shared Energy of LayerCluster vs #eta per best Calo Particle for layer " + istr2,
0775                           nintEta_,
0776                           minEta_,
0777                           maxEta_,
0778                           minSharedEneFrac_,
0779                           maxSharedEneFrac_);
0780     histograms.h_sharedenergy_layercl2caloparticle_vs_phi_perlayer[ilayer] =
0781         ibook.bookProfile("SharedEnergy_layercl2caloparticle_vs_phi_perlayer" + istr1,
0782                           "Shared Energy of LayerCluster vs #phi per best Calo Particle for layer " + istr2,
0783                           nintPhi_,
0784                           minPhi_,
0785                           maxPhi_,
0786                           minSharedEneFrac_,
0787                           maxSharedEneFrac_);
0788     histograms.h_num_caloparticle_eta_perlayer[ilayer] =
0789         ibook.book1D("Num_CaloParticle_Eta_perlayer" + istr1,
0790                      "Num CaloParticle Eta per Layer Cluster for layer " + istr2,
0791                      nintEta_,
0792                      minEta_,
0793                      maxEta_);
0794     histograms.h_numDup_caloparticle_eta_perlayer[ilayer] =
0795         ibook.book1D("NumDup_CaloParticle_Eta_perlayer" + istr1,
0796                      "Num Duplicate CaloParticle Eta per Layer Cluster for layer " + istr2,
0797                      nintEta_,
0798                      minEta_,
0799                      maxEta_);
0800     histograms.h_denom_caloparticle_eta_perlayer[ilayer] =
0801         ibook.book1D("Denom_CaloParticle_Eta_perlayer" + istr1,
0802                      "Denom CaloParticle Eta per Layer Cluster for layer " + istr2,
0803                      nintEta_,
0804                      minEta_,
0805                      maxEta_);
0806     histograms.h_num_caloparticle_phi_perlayer[ilayer] =
0807         ibook.book1D("Num_CaloParticle_Phi_perlayer" + istr1,
0808                      "Num CaloParticle Phi per Layer Cluster for layer " + istr2,
0809                      nintPhi_,
0810                      minPhi_,
0811                      maxPhi_);
0812     histograms.h_numDup_caloparticle_phi_perlayer[ilayer] =
0813         ibook.book1D("NumDup_CaloParticle_Phi_perlayer" + istr1,
0814                      "Num Duplicate CaloParticle Phi per Layer Cluster for layer " + istr2,
0815                      nintPhi_,
0816                      minPhi_,
0817                      maxPhi_);
0818     histograms.h_denom_caloparticle_phi_perlayer[ilayer] =
0819         ibook.book1D("Denom_CaloParticle_Phi_perlayer" + istr1,
0820                      "Denom CaloParticle Phi per Layer Cluster for layer " + istr2,
0821                      nintPhi_,
0822                      minPhi_,
0823                      maxPhi_);
0824     histograms.h_num_layercl_eta_perlayer[ilayer] =
0825         ibook.book1D("Num_LayerCluster_Eta_perlayer" + istr1,
0826                      "Num LayerCluster Eta per Layer Cluster for layer " + istr2,
0827                      nintEta_,
0828                      minEta_,
0829                      maxEta_);
0830     histograms.h_numMerge_layercl_eta_perlayer[ilayer] =
0831         ibook.book1D("NumMerge_LayerCluster_Eta_perlayer" + istr1,
0832                      "Num Merge LayerCluster Eta per Layer Cluster for layer " + istr2,
0833                      nintEta_,
0834                      minEta_,
0835                      maxEta_);
0836     histograms.h_denom_layercl_eta_perlayer[ilayer] =
0837         ibook.book1D("Denom_LayerCluster_Eta_perlayer" + istr1,
0838                      "Denom LayerCluster Eta per Layer Cluster for layer " + istr2,
0839                      nintEta_,
0840                      minEta_,
0841                      maxEta_);
0842     histograms.h_num_layercl_phi_perlayer[ilayer] =
0843         ibook.book1D("Num_LayerCluster_Phi_perlayer" + istr1,
0844                      "Num LayerCluster Phi per Layer Cluster for layer " + istr2,
0845                      nintPhi_,
0846                      minPhi_,
0847                      maxPhi_);
0848     histograms.h_numMerge_layercl_phi_perlayer[ilayer] =
0849         ibook.book1D("NumMerge_LayerCluster_Phi_perlayer" + istr1,
0850                      "Num Merge LayerCluster Phi per Layer Cluster for layer " + istr2,
0851                      nintPhi_,
0852                      minPhi_,
0853                      maxPhi_);
0854     histograms.h_denom_layercl_phi_perlayer[ilayer] =
0855         ibook.book1D("Denom_LayerCluster_Phi_perlayer" + istr1,
0856                      "Denom LayerCluster Phi per Layer Cluster for layer " + istr2,
0857                      nintPhi_,
0858                      minPhi_,
0859                      maxPhi_);
0860   }
0861   //---------------------------------------------------------------------------------------------------------------------------
0862 }
0863 
0864 void HGVHistoProducerAlgo::bookClusterHistos_CellLevel(DQMStore::IBooker& ibook,
0865                                                        Histograms& histograms,
0866                                                        unsigned int layers,
0867                                                        std::vector<int> thicknesses) {
0868   //----------------------------------------------------------------------------------------------------------------------------
0869   for (unsigned ilayer = 0; ilayer < 2 * layers; ++ilayer) {
0870     auto istr1 = std::to_string(ilayer);
0871     while (istr1.size() < 2) {
0872       istr1.insert(0, "0");
0873     }
0874     // Make a mapping to the regural layer naming plus z- or z+ for convenience
0875     std::string istr2 = "";
0876     // first with the -z endcap
0877     if (ilayer < layers) {
0878       istr2 = std::to_string(ilayer + 1) + " in z-";
0879     } else {  // then for the +z
0880       istr2 = std::to_string(ilayer - (layers - 1)) + " in z+";
0881     }
0882     histograms.h_cellAssociation_perlayer[ilayer] =
0883         ibook.book1D("cellAssociation_perlayer" + istr1, "Cell Association for layer " + istr2, 5, -4., 1.);
0884     histograms.h_cellAssociation_perlayer[ilayer]->setBinLabel(2, "TN(purity)");
0885     histograms.h_cellAssociation_perlayer[ilayer]->setBinLabel(3, "FN(ineff.)");
0886     histograms.h_cellAssociation_perlayer[ilayer]->setBinLabel(4, "FP(fake)");
0887     histograms.h_cellAssociation_perlayer[ilayer]->setBinLabel(5, "TP(eff.)");
0888   }
0889   //----------------------------------------------------------------------------------------------------------------------------
0890   for (std::vector<int>::iterator it = thicknesses.begin(); it != thicknesses.end(); ++it) {
0891     auto istr = std::to_string(*it);
0892     histograms.h_cellsenedens_perthick[(*it)] = ibook.book1D("cellsenedens_thick_" + istr,
0893                                                              "energy density of cluster cells for thickness " + istr,
0894                                                              nintCellsEneDensperthick_,
0895                                                              minCellsEneDensperthick_,
0896                                                              maxCellsEneDensperthick_);
0897   }
0898   //----------------------------------------------------------------------------------------------------------------------------
0899   //Not all combination exists but should keep them all for cross checking reason.
0900   for (std::vector<int>::iterator it = thicknesses.begin(); it != thicknesses.end(); ++it) {
0901     for (unsigned ilayer = 0; ilayer < 2 * layers; ++ilayer) {
0902       auto istr1 = std::to_string(*it);
0903       auto istr2 = std::to_string(ilayer);
0904       while (istr2.size() < 2)
0905         istr2.insert(0, "0");
0906       auto istr = istr1 + "_" + istr2;
0907       // Make a mapping to the regural layer naming plus z- or z+ for convenience
0908       std::string istr3 = "";
0909       // first with the -z endcap
0910       if (ilayer < layers) {
0911         istr3 = std::to_string(ilayer + 1) + " in z- ";
0912       } else {  // then for the +z
0913         istr3 = std::to_string(ilayer - (layers - 1)) + " in z+ ";
0914       }
0915       //---
0916       histograms.h_cellsnum_perthickperlayer[istr] =
0917           ibook.book1D("cellsnum_perthick_perlayer_" + istr,
0918                        "total number of cells for layer " + istr3 + " for thickness " + istr1,
0919                        nintTotNcellsperthickperlayer_,
0920                        minTotNcellsperthickperlayer_,
0921                        maxTotNcellsperthickperlayer_);
0922       //---
0923       histograms.h_distancetoseedcell_perthickperlayer[istr] =
0924           ibook.book1D("distancetoseedcell_perthickperlayer_" + istr,
0925                        "distance of cluster cells to seed cell for layer " + istr3 + " for thickness " + istr1,
0926                        nintDisToSeedperthickperlayer_,
0927                        minDisToSeedperthickperlayer_,
0928                        maxDisToSeedperthickperlayer_);
0929       //---
0930       histograms.h_distancetoseedcell_perthickperlayer_eneweighted[istr] = ibook.book1D(
0931           "distancetoseedcell_perthickperlayer_eneweighted_" + istr,
0932           "energy weighted distance of cluster cells to seed cell for layer " + istr3 + " for thickness " + istr1,
0933           nintDisToSeedperthickperlayerenewei_,
0934           minDisToSeedperthickperlayerenewei_,
0935           maxDisToSeedperthickperlayerenewei_);
0936       //---
0937       histograms.h_distancetomaxcell_perthickperlayer[istr] =
0938           ibook.book1D("distancetomaxcell_perthickperlayer_" + istr,
0939                        "distance of cluster cells to max cell for layer " + istr3 + " for thickness " + istr1,
0940                        nintDisToMaxperthickperlayer_,
0941                        minDisToMaxperthickperlayer_,
0942                        maxDisToMaxperthickperlayer_);
0943       //---
0944       histograms.h_distancetomaxcell_perthickperlayer_eneweighted[istr] = ibook.book1D(
0945           "distancetomaxcell_perthickperlayer_eneweighted_" + istr,
0946           "energy weighted distance of cluster cells to max cell for layer " + istr3 + " for thickness " + istr1,
0947           nintDisToMaxperthickperlayerenewei_,
0948           minDisToMaxperthickperlayerenewei_,
0949           maxDisToMaxperthickperlayerenewei_);
0950       //---
0951       histograms.h_distancebetseedandmaxcell_perthickperlayer[istr] =
0952           ibook.book1D("distancebetseedandmaxcell_perthickperlayer_" + istr,
0953                        "distance of seed cell to max cell for layer " + istr3 + " for thickness " + istr1,
0954                        nintDisSeedToMaxperthickperlayer_,
0955                        minDisSeedToMaxperthickperlayer_,
0956                        maxDisSeedToMaxperthickperlayer_);
0957       //---
0958       histograms.h_distancebetseedandmaxcellvsclusterenergy_perthickperlayer[istr] = ibook.book2D(
0959           "distancebetseedandmaxcellvsclusterenergy_perthickperlayer_" + istr,
0960           "distance of seed cell to max cell vs cluster energy for layer " + istr3 + " for thickness " + istr1,
0961           nintDisSeedToMaxperthickperlayer_,
0962           minDisSeedToMaxperthickperlayer_,
0963           maxDisSeedToMaxperthickperlayer_,
0964           nintClEneperthickperlayer_,
0965           minClEneperthickperlayer_,
0966           maxClEneperthickperlayer_);
0967     }
0968   }
0969 }
0970 //----------------------------------------------------------------------------------------------------------------------------
0971 
0972 void HGVHistoProducerAlgo::bookTracksterHistos(DQMStore::IBooker& ibook, Histograms& histograms, unsigned int layers) {
0973   std::unordered_map<int, dqm::reco::MonitorElement*> clusternum_in_trackster_perlayer;
0974   clusternum_in_trackster_perlayer.clear();
0975 
0976   for (unsigned ilayer = 0; ilayer < 2 * layers; ++ilayer) {
0977     auto istr1 = std::to_string(ilayer);
0978     while (istr1.size() < 2) {
0979       istr1.insert(0, "0");
0980     }
0981     // Make a mapping to the regural layer naming plus z- or z+ for convenience
0982     std::string istr2 = "";
0983     // first with the -z endcap
0984     if (ilayer < layers) {
0985       istr2 = std::to_string(ilayer + 1) + " in z-";
0986     } else {  // then for the +z
0987       istr2 = std::to_string(ilayer - (layers - 1)) + " in z+";
0988     }
0989 
0990     clusternum_in_trackster_perlayer[ilayer] = ibook.book1D("clusternum_in_trackster_perlayer" + istr1,
0991                                                             "Number of layer clusters in Trackster for layer " + istr2,
0992                                                             nintTotNClsinTSTsperlayer_,
0993                                                             minTotNClsinTSTsperlayer_,
0994                                                             maxTotNClsinTSTsperlayer_);
0995   }
0996 
0997   histograms.h_clusternum_in_trackster_perlayer.push_back(std::move(clusternum_in_trackster_perlayer));
0998 
0999   histograms.h_tracksternum.push_back(ibook.book1D(
1000       "tottracksternum", "total number of Tracksters;# of Tracksters", nintTotNTSTs_, minTotNTSTs_, maxTotNTSTs_));
1001 
1002   histograms.h_conttracksternum.push_back(ibook.book1D(
1003       "conttracksternum", "number of Tracksters with 3 contiguous layers", nintTotNTSTs_, minTotNTSTs_, maxTotNTSTs_));
1004 
1005   histograms.h_nonconttracksternum.push_back(ibook.book1D("nonconttracksternum",
1006                                                           "number of Tracksters without 3 contiguous layers",
1007                                                           nintTotNTSTs_,
1008                                                           minTotNTSTs_,
1009                                                           maxTotNTSTs_));
1010 
1011   histograms.h_clusternum_in_trackster.push_back(
1012       ibook.book1D("clusternum_in_trackster",
1013                    "total number of layer clusters in Trackster;# of LayerClusters",
1014                    nintTotNClsinTSTs_,
1015                    minTotNClsinTSTs_,
1016                    maxTotNClsinTSTs_));
1017 
1018   histograms.h_clusternum_in_trackster_vs_layer.push_back(ibook.bookProfile(
1019       "clusternum_in_trackster_vs_layer",
1020       "Profile of 2d layer clusters in Trackster vs layer number;layer number;<2D LayerClusters in Trackster>",
1021       2 * layers,
1022       0.,
1023       2. * layers,
1024       minTotNClsinTSTsperlayer_,
1025       maxTotNClsinTSTsperlayer_));
1026 
1027   histograms.h_multiplicityOfLCinTST.push_back(
1028       ibook.book2D("multiplicityOfLCinTST",
1029                    "Multiplicity vs Layer cluster size in Tracksters;LayerCluster multiplicity in Tracksters;Cluster "
1030                    "size (n_{hits})",
1031                    nintMplofLCs_,
1032                    minMplofLCs_,
1033                    maxMplofLCs_,
1034                    nintSizeCLsinTSTs_,
1035                    minSizeCLsinTSTs_,
1036                    maxSizeCLsinTSTs_));
1037 
1038   histograms.h_multiplicity_numberOfEventsHistogram.push_back(ibook.book1D("multiplicity_numberOfEventsHistogram",
1039                                                                            "multiplicity numberOfEventsHistogram",
1040                                                                            nintMplofLCs_,
1041                                                                            minMplofLCs_,
1042                                                                            maxMplofLCs_));
1043 
1044   histograms.h_multiplicity_zminus_numberOfEventsHistogram.push_back(
1045       ibook.book1D("multiplicity_zminus_numberOfEventsHistogram",
1046                    "multiplicity numberOfEventsHistogram in z-",
1047                    nintMplofLCs_,
1048                    minMplofLCs_,
1049                    maxMplofLCs_));
1050 
1051   histograms.h_multiplicity_zplus_numberOfEventsHistogram.push_back(
1052       ibook.book1D("multiplicity_zplus_numberOfEventsHistogram",
1053                    "multiplicity numberOfEventsHistogram in z+",
1054                    nintMplofLCs_,
1055                    minMplofLCs_,
1056                    maxMplofLCs_));
1057 
1058   histograms.h_multiplicityOfLCinTST_vs_layercluster_zminus.push_back(
1059       ibook.book2D("multiplicityOfLCinTST_vs_layercluster_zminus",
1060                    "Multiplicity vs Layer number in z-;LayerCluster multiplicity in Tracksters;layer number",
1061                    nintMplofLCs_,
1062                    minMplofLCs_,
1063                    maxMplofLCs_,
1064                    layers,
1065                    0.,
1066                    (float)layers));
1067 
1068   histograms.h_multiplicityOfLCinTST_vs_layercluster_zplus.push_back(
1069       ibook.book2D("multiplicityOfLCinTST_vs_layercluster_zplus",
1070                    "Multiplicity vs Layer number in z+;LayerCluster multiplicity in Tracksters;layer number",
1071                    nintMplofLCs_,
1072                    minMplofLCs_,
1073                    maxMplofLCs_,
1074                    layers,
1075                    0.,
1076                    (float)layers));
1077 
1078   histograms.h_multiplicityOfLCinTST_vs_layerclusterenergy.push_back(
1079       ibook.book2D("multiplicityOfLCinTST_vs_layerclusterenergy",
1080                    "Multiplicity vs Layer cluster energy;LayerCluster multiplicity in Tracksters;Cluster energy [GeV]",
1081                    nintMplofLCs_,
1082                    minMplofLCs_,
1083                    maxMplofLCs_,
1084                    nintClEnepermultiplicity_,
1085                    minClEnepermultiplicity_,
1086                    maxClEnepermultiplicity_));
1087 
1088   histograms.h_trackster_pt.push_back(
1089       ibook.book1D("trackster_pt", "Pt of the Trackster;Trackster p_{T} [GeV]", nintPt_, minPt_, maxPt_));
1090   histograms.h_trackster_eta.push_back(
1091       ibook.book1D("trackster_eta", "Eta of the Trackster;Trackster #eta", nintEta_, minEta_, maxEta_));
1092   histograms.h_trackster_phi.push_back(
1093       ibook.book1D("trackster_phi", "Phi of the Trackster;Trackster #phi", nintPhi_, minPhi_, maxPhi_));
1094   histograms.h_trackster_energy.push_back(
1095       ibook.book1D("trackster_energy", "Energy of the Trackster;Trackster energy [GeV]", nintEne_, minEne_, maxEne_));
1096   histograms.h_trackster_x.push_back(
1097       ibook.book1D("trackster_x", "X position of the Trackster;Trackster x", nintX_, minX_, maxX_));
1098   histograms.h_trackster_y.push_back(
1099       ibook.book1D("trackster_y", "Y position of the Trackster;Trackster y", nintY_, minY_, maxY_));
1100   histograms.h_trackster_z.push_back(
1101       ibook.book1D("trackster_z", "Z position of the Trackster;Trackster z", nintZ_, minZ_, maxZ_));
1102   histograms.h_trackster_firstlayer.push_back(ibook.book1D(
1103       "trackster_firstlayer", "First layer of the Trackster;Trackster First Layer", 2 * layers, 0., (float)2 * layers));
1104   histograms.h_trackster_lastlayer.push_back(ibook.book1D(
1105       "trackster_lastlayer", "Last layer of the Trackster;Trackster Last Layer", 2 * layers, 0., (float)2 * layers));
1106   histograms.h_trackster_layersnum.push_back(
1107       ibook.book1D("trackster_layersnum",
1108                    "Number of layers of the Trackster;Trackster Number of Layers",
1109                    2 * layers,
1110                    0.,
1111                    (float)2 * layers));
1112 }
1113 
1114 void HGVHistoProducerAlgo::bookTracksterSTSHistos(DQMStore::IBooker& ibook,
1115                                                   Histograms& histograms,
1116                                                   const validationType valType) {
1117   const string ref[] = {"caloparticle", "simtrackster", "simtrackster_fromCP"};
1118   const string refT[] = {"CaloParticle", "SimTrackster", "SimTrackster_fromCP"};
1119   // Must be in sync with labels in PostProcessorHGCAL_cfi.py
1120   const string val[] = {"_Link", "_PR"};
1121   //const string val[] = {"_CP", "_PR", "_Link"};
1122   const string rtos = ";score Reco-to-Sim";
1123   const string stor = ";score Sim-to-Reco";
1124   const string shREnFr = ";shared Reco energy fraction";
1125   const string shSEnFr = ";shared Sim energy fraction";
1126 
1127   histograms.h_score_trackster2caloparticle[valType].push_back(
1128       ibook.book1D("Score_trackster2" + ref[valType],
1129                    "Score of Trackster per " + refT[valType] + rtos,
1130                    nintScore_,
1131                    minScore_,
1132                    maxScore_));
1133   histograms.h_score_trackster2bestCaloparticle[valType].push_back(
1134       ibook.book1D("ScoreFake_trackster2" + ref[valType],
1135                    "Score of Trackster per best " + refT[valType] + rtos,
1136                    nintScore_,
1137                    minScore_,
1138                    maxScore_));
1139   histograms.h_score_trackster2bestCaloparticle2[valType].push_back(
1140       ibook.book1D("ScoreMerge_trackster2" + ref[valType],
1141                    "Score of Trackster per 2^{nd} best " + refT[valType] + rtos,
1142                    nintScore_,
1143                    minScore_,
1144                    maxScore_));
1145   histograms.h_score_caloparticle2trackster[valType].push_back(
1146       ibook.book1D("Score_" + ref[valType] + "2trackster",
1147                    "Score of " + refT[valType] + " per Trackster" + stor,
1148                    nintScore_,
1149                    minScore_,
1150                    maxScore_));
1151   histograms.h_scorePur_caloparticle2trackster[valType].push_back(
1152       ibook.book1D("ScorePur_" + ref[valType] + "2trackster",
1153                    "Score of " + refT[valType] + " per best Trackster" + stor,
1154                    nintScore_,
1155                    minScore_,
1156                    maxScore_));
1157   histograms.h_scoreDupl_caloparticle2trackster[valType].push_back(
1158       ibook.book1D("ScoreDupl_" + ref[valType] + "2trackster",
1159                    "Score of " + refT[valType] + " per 2^{nd} best Trackster" + stor,
1160                    nintScore_,
1161                    minScore_,
1162                    maxScore_));
1163   histograms.h_energy_vs_score_trackster2caloparticle[valType].push_back(
1164       ibook.book2D("Energy_vs_Score_trackster2" + refT[valType],
1165                    "Energy vs Score of Trackster per " + refT[valType] + rtos + shREnFr,
1166                    nintScore_,
1167                    minScore_,
1168                    maxScore_,
1169                    nintSharedEneFrac_,
1170                    minTSTSharedEneFrac_,
1171                    maxTSTSharedEneFrac_));
1172   histograms.h_energy_vs_score_trackster2bestCaloparticle[valType].push_back(
1173       ibook.book2D("Energy_vs_Score_trackster2best" + refT[valType],
1174                    "Energy vs Score of Trackster per best " + refT[valType] + rtos + shREnFr,
1175                    nintScore_,
1176                    minScore_,
1177                    maxScore_,
1178                    nintSharedEneFrac_,
1179                    minTSTSharedEneFrac_,
1180                    maxTSTSharedEneFrac_));
1181   histograms.h_energy_vs_score_trackster2bestCaloparticle2[valType].push_back(
1182       ibook.book2D("Energy_vs_Score_trackster2secBest" + refT[valType],
1183                    "Energy vs Score of Trackster per 2^{nd} best " + refT[valType] + rtos + shREnFr,
1184                    nintScore_,
1185                    minScore_,
1186                    maxScore_,
1187                    nintSharedEneFrac_,
1188                    minTSTSharedEneFrac_,
1189                    maxTSTSharedEneFrac_));
1190   histograms.h_energy_vs_score_caloparticle2trackster[valType].push_back(
1191       ibook.book2D("Energy_vs_Score_" + ref[valType] + "2Trackster",
1192                    "Energy vs Score of " + refT[valType] + " per Trackster" + stor + shSEnFr,
1193                    nintScore_,
1194                    minScore_,
1195                    maxScore_,
1196                    nintSharedEneFrac_,
1197                    minTSTSharedEneFrac_,
1198                    maxTSTSharedEneFrac_));
1199   histograms.h_energy_vs_score_caloparticle2bestTrackster[valType].push_back(
1200       ibook.book2D("Energy_vs_Score_" + ref[valType] + "2bestTrackster",
1201                    "Energy vs Score of " + refT[valType] + " per best Trackster" + stor + shSEnFr,
1202                    nintScore_,
1203                    minScore_,
1204                    maxScore_,
1205                    nintSharedEneFrac_,
1206                    minTSTSharedEneFrac_,
1207                    maxTSTSharedEneFrac_));
1208   histograms.h_energy_vs_score_caloparticle2bestTrackster2[valType].push_back(
1209       ibook.book2D("Energy_vs_Score_" + ref[valType] + "2secBestTrackster",
1210                    "Energy vs Score of " + refT[valType] + " per 2^{nd} best Trackster" + stor + shSEnFr,
1211                    nintScore_,
1212                    minScore_,
1213                    maxScore_,
1214                    nintSharedEneFrac_,
1215                    minTSTSharedEneFrac_,
1216                    maxTSTSharedEneFrac_));
1217 
1218   // Back to all Tracksters
1219   // eta
1220   histograms.h_num_trackster_eta[valType].push_back(ibook.book1D(
1221       "Num_Trackster_Eta" + val[valType], "Num Trackster Eta per Trackster;#eta", nintEta_, minEta_, maxEta_));
1222   histograms.h_numMerge_trackster_eta[valType].push_back(ibook.book1D("NumMerge_Trackster_Eta" + val[valType],
1223                                                                       "Num Merge Trackster Eta per Trackster;#eta",
1224                                                                       nintEta_,
1225                                                                       minEta_,
1226                                                                       maxEta_));
1227   histograms.h_denom_trackster_eta[valType].push_back(ibook.book1D(
1228       "Denom_Trackster_Eta" + val[valType], "Denom Trackster Eta per Trackster;#eta", nintEta_, minEta_, maxEta_));
1229   // phi
1230   histograms.h_num_trackster_phi[valType].push_back(ibook.book1D(
1231       "Num_Trackster_Phi" + val[valType], "Num Trackster Phi per Trackster;#phi", nintPhi_, minPhi_, maxPhi_));
1232   histograms.h_numMerge_trackster_phi[valType].push_back(ibook.book1D("NumMerge_Trackster_Phi" + val[valType],
1233                                                                       "Num Merge Trackster Phi per Trackster;#phi",
1234                                                                       nintPhi_,
1235                                                                       minPhi_,
1236                                                                       maxPhi_));
1237   histograms.h_denom_trackster_phi[valType].push_back(ibook.book1D(
1238       "Denom_Trackster_Phi" + val[valType], "Denom Trackster Phi per Trackster;#phi", nintPhi_, minPhi_, maxPhi_));
1239   // energy
1240   histograms.h_num_trackster_en[valType].push_back(ibook.book1D("Num_Trackster_Energy" + val[valType],
1241                                                                 "Num Trackster Energy per Trackster;energy [GeV]",
1242                                                                 nintEne_,
1243                                                                 minEne_,
1244                                                                 maxEne_));
1245   histograms.h_numMerge_trackster_en[valType].push_back(
1246       ibook.book1D("NumMerge_Trackster_Energy" + val[valType],
1247                    "Num Merge Trackster Energy per Trackster;energy [GeV]",
1248                    nintEne_,
1249                    minEne_,
1250                    maxEne_));
1251   histograms.h_denom_trackster_en[valType].push_back(ibook.book1D("Denom_Trackster_Energy" + val[valType],
1252                                                                   "Denom Trackster Energy per Trackster;energy [GeV]",
1253                                                                   nintEne_,
1254                                                                   minEne_,
1255                                                                   maxEne_));
1256   // pT
1257   histograms.h_num_trackster_pt[valType].push_back(ibook.book1D(
1258       "Num_Trackster_Pt" + val[valType], "Num Trackster p_{T} per Trackster;p_{T} [GeV]", nintPt_, minPt_, maxPt_));
1259   histograms.h_numMerge_trackster_pt[valType].push_back(
1260       ibook.book1D("NumMerge_Trackster_Pt" + val[valType],
1261                    "Num Merge Trackster p_{T} per Trackster;p_{T} [GeV]",
1262                    nintPt_,
1263                    minPt_,
1264                    maxPt_));
1265   histograms.h_denom_trackster_pt[valType].push_back(ibook.book1D(
1266       "Denom_Trackster_Pt" + val[valType], "Denom Trackster p_{T} per Trackster;p_{T} [GeV]", nintPt_, minPt_, maxPt_));
1267 
1268   histograms.h_sharedenergy_trackster2caloparticle[valType].push_back(
1269       ibook.book1D("SharedEnergy_trackster2" + ref[valType],
1270                    "Shared Energy of Trackster per " + refT[valType] + shREnFr,
1271                    nintSharedEneFrac_,
1272                    minTSTSharedEneFrac_,
1273                    maxTSTSharedEneFrac_));
1274   histograms.h_sharedenergy_trackster2bestCaloparticle[valType].push_back(
1275       ibook.book1D("SharedEnergy_trackster2" + ref[valType] + "_assoc",
1276                    "Shared Energy of Trackster per best " + refT[valType] + shREnFr,
1277                    nintSharedEneFrac_,
1278                    minTSTSharedEneFrac_,
1279                    maxTSTSharedEneFrac_));
1280   histograms.h_sharedenergy_trackster2bestCaloparticle_vs_eta[valType].push_back(
1281       ibook.bookProfile("SharedEnergy_trackster2" + ref[valType] + "_assoc_vs_eta",
1282                         "Shared Energy of Trackster vs #eta per best " + refT[valType] + ";Trackster #eta" + shREnFr,
1283                         nintEta_,
1284                         minEta_,
1285                         maxEta_,
1286                         minTSTSharedEneFrac_,
1287                         maxTSTSharedEneFrac_));
1288   histograms.h_sharedenergy_trackster2bestCaloparticle_vs_phi[valType].push_back(
1289       ibook.bookProfile("SharedEnergy_trackster2" + ref[valType] + "_assoc_vs_phi",
1290                         "Shared Energy of Trackster vs #phi per best " + refT[valType] + ";Trackster #phi" + shREnFr,
1291                         nintPhi_,
1292                         minPhi_,
1293                         maxPhi_,
1294                         minTSTSharedEneFrac_,
1295                         maxTSTSharedEneFrac_));
1296   histograms.h_sharedenergy_trackster2bestCaloparticle2[valType].push_back(
1297       ibook.book1D("SharedEnergy_trackster2" + ref[valType] + "_assoc2",
1298                    "Shared Energy of Trackster per 2^{nd} best " + refT[valType] + shREnFr,
1299                    nintSharedEneFrac_,
1300                    minTSTSharedEneFrac_,
1301                    maxTSTSharedEneFrac_));
1302 
1303   histograms.h_sharedenergy_caloparticle2trackster[valType].push_back(
1304       ibook.book1D("SharedEnergy_" + ref[valType] + "2trackster",
1305                    "Shared Energy of " + refT[valType] + " per Trackster" + shSEnFr,
1306                    nintSharedEneFrac_,
1307                    minTSTSharedEneFrac_,
1308                    maxTSTSharedEneFrac_));
1309   histograms.h_sharedenergy_caloparticle2trackster_assoc[valType].push_back(
1310       ibook.book1D("SharedEnergy_" + ref[valType] + "2trackster_assoc",
1311                    "Shared Energy of " + refT[valType] + " per best Trackster" + shSEnFr,
1312                    nintSharedEneFrac_,
1313                    minTSTSharedEneFrac_,
1314                    maxTSTSharedEneFrac_));
1315   histograms.h_sharedenergy_caloparticle2trackster_assoc_vs_eta[valType].push_back(ibook.bookProfile(
1316       "SharedEnergy_" + ref[valType] + "2trackster_assoc_vs_eta",
1317       "Shared Energy of " + refT[valType] + " vs #eta per best Trackster;" + refT[valType] + " #eta" + shSEnFr,
1318       nintEta_,
1319       minEta_,
1320       maxEta_,
1321       minTSTSharedEneFrac_,
1322       maxTSTSharedEneFrac_));
1323   histograms.h_sharedenergy_caloparticle2trackster_assoc_vs_phi[valType].push_back(ibook.bookProfile(
1324       "SharedEnergy_" + ref[valType] + "2trackster_assoc_vs_phi",
1325       "Shared Energy of " + refT[valType] + " vs #phi per best Trackster;" + refT[valType] + " #phi" + shSEnFr,
1326       nintPhi_,
1327       minPhi_,
1328       maxPhi_,
1329       minTSTSharedEneFrac_,
1330       maxTSTSharedEneFrac_));
1331   histograms.h_sharedenergy_caloparticle2trackster_assoc2[valType].push_back(
1332       ibook.book1D("SharedEnergy_" + ref[valType] + "2trackster_assoc2",
1333                    "Shared Energy of " + refT[valType] + " per 2^{nd} best Trackster;" + shSEnFr,
1334                    nintSharedEneFrac_,
1335                    minTSTSharedEneFrac_,
1336                    maxTSTSharedEneFrac_));
1337 
1338   // eta
1339   histograms.h_numEff_caloparticle_eta[valType].push_back(
1340       ibook.book1D("NumEff_" + refT[valType] + "_Eta",
1341                    "Num Efficiency " + refT[valType] + " Eta per Trackster;#eta",
1342                    nintEta_,
1343                    minEta_,
1344                    maxEta_));
1345   histograms.h_num_caloparticle_eta[valType].push_back(
1346       ibook.book1D("Num_" + refT[valType] + "_Eta",
1347                    "Num Purity " + refT[valType] + " Eta per Trackster;#eta",
1348                    nintEta_,
1349                    minEta_,
1350                    maxEta_));
1351   histograms.h_numDup_trackster_eta[valType].push_back(ibook.book1D(
1352       "NumDup_Trackster_Eta" + val[valType], "Num Duplicate Trackster vs Eta;#eta", nintEta_, minEta_, maxEta_));
1353   histograms.h_denom_caloparticle_eta[valType].push_back(
1354       ibook.book1D("Denom_" + refT[valType] + "_Eta",
1355                    "Denom " + refT[valType] + " Eta per Trackster;#eta",
1356                    nintEta_,
1357                    minEta_,
1358                    maxEta_));
1359   // phi
1360   histograms.h_numEff_caloparticle_phi[valType].push_back(
1361       ibook.book1D("NumEff_" + refT[valType] + "_Phi",
1362                    "Num Efficiency " + refT[valType] + " Phi per Trackster;#phi",
1363                    nintPhi_,
1364                    minPhi_,
1365                    maxPhi_));
1366   histograms.h_num_caloparticle_phi[valType].push_back(
1367       ibook.book1D("Num_" + refT[valType] + "_Phi",
1368                    "Num Purity " + refT[valType] + " Phi per Trackster;#phi",
1369                    nintPhi_,
1370                    minPhi_,
1371                    maxPhi_));
1372   histograms.h_numDup_trackster_phi[valType].push_back(ibook.book1D(
1373       "NumDup_Trackster_Phi" + val[valType], "Num Duplicate Trackster vs Phi;#phi", nintPhi_, minPhi_, maxPhi_));
1374   histograms.h_denom_caloparticle_phi[valType].push_back(
1375       ibook.book1D("Denom_" + refT[valType] + "_Phi",
1376                    "Denom " + refT[valType] + " Phi per Trackster;#phi",
1377                    nintPhi_,
1378                    minPhi_,
1379                    maxPhi_));
1380   // energy
1381   histograms.h_numEff_caloparticle_en[valType].push_back(
1382       ibook.book1D("NumEff_" + refT[valType] + "_Energy",
1383                    "Num Efficiency " + refT[valType] + " Energy per Trackster;energy [GeV]",
1384                    nintEne_,
1385                    minEne_,
1386                    maxEne_));
1387   histograms.h_num_caloparticle_en[valType].push_back(
1388       ibook.book1D("Num_" + refT[valType] + "_Energy",
1389                    "Num Purity " + refT[valType] + " Energy per Trackster;energy [GeV]",
1390                    nintEne_,
1391                    minEne_,
1392                    maxEne_));
1393   histograms.h_numDup_trackster_en[valType].push_back(ibook.book1D("NumDup_Trackster_Energy" + val[valType],
1394                                                                    "Num Duplicate Trackster vs Energy;energy [GeV]",
1395                                                                    nintEne_,
1396                                                                    minEne_,
1397                                                                    maxEne_));
1398   histograms.h_denom_caloparticle_en[valType].push_back(
1399       ibook.book1D("Denom_" + refT[valType] + "_Energy",
1400                    "Denom " + refT[valType] + " Energy per Trackster;energy [GeV]",
1401                    nintEne_,
1402                    minEne_,
1403                    maxEne_));
1404   // pT
1405   histograms.h_numEff_caloparticle_pt[valType].push_back(
1406       ibook.book1D("NumEff_" + refT[valType] + "_Pt",
1407                    "Num Efficiency " + refT[valType] + " p_{T} per Trackster;p_{T} [GeV]",
1408                    nintPt_,
1409                    minPt_,
1410                    maxPt_));
1411   histograms.h_num_caloparticle_pt[valType].push_back(
1412       ibook.book1D("Num_" + refT[valType] + "_Pt",
1413                    "Num Purity " + refT[valType] + " p_{T} per Trackster;p_{T} [GeV]",
1414                    nintPt_,
1415                    minPt_,
1416                    maxPt_));
1417   histograms.h_numDup_trackster_pt[valType].push_back(ibook.book1D(
1418       "NumDup_Trackster_Pt" + val[valType], "Num Duplicate Trackster vs p_{T};p_{T} [GeV]", nintPt_, minPt_, maxPt_));
1419   histograms.h_denom_caloparticle_pt[valType].push_back(
1420       ibook.book1D("Denom_" + refT[valType] + "_Pt",
1421                    "Denom " + refT[valType] + " p_{T} per Trackster;p_{T} [GeV]",
1422                    nintPt_,
1423                    minPt_,
1424                    maxPt_));
1425 }
1426 
1427 void HGVHistoProducerAlgo::fill_info_histos(const Histograms& histograms, unsigned int layers) const {
1428   // Save some info straight from geometry to avoid mistakes from updates
1429   //----------- TODO ----------------------------------------------------------
1430   // For now values returned for 'lastLayerFHzp': '104', 'lastLayerFHzm': '52' are not the one expected.
1431   // Will come back to this when there will be info in CMSSW to put in DQM file.
1432   histograms.lastLayerEEzm->Fill(recHitTools_->lastLayerEE());
1433   histograms.lastLayerFHzm->Fill(recHitTools_->lastLayerFH());
1434   histograms.maxlayerzm->Fill(layers);
1435   histograms.lastLayerEEzp->Fill(recHitTools_->lastLayerEE() + layers);
1436   histograms.lastLayerFHzp->Fill(recHitTools_->lastLayerFH() + layers);
1437   histograms.maxlayerzp->Fill(layers + layers);
1438 }
1439 
1440 void HGVHistoProducerAlgo::fill_caloparticle_histos(const Histograms& histograms,
1441                                                     int pdgid,
1442                                                     const CaloParticle& caloParticle,
1443                                                     std::vector<SimVertex> const& simVertices,
1444                                                     unsigned int layers,
1445                                                     std::unordered_map<DetId, const HGCRecHit*> const& hitMap) const {
1446   const auto eta = getEta(caloParticle.eta());
1447   if (histograms.h_caloparticle_eta.count(pdgid)) {
1448     histograms.h_caloparticle_eta.at(pdgid)->Fill(eta);
1449   }
1450   if (histograms.h_caloparticle_eta_Zorigin.count(pdgid)) {
1451     histograms.h_caloparticle_eta_Zorigin.at(pdgid)->Fill(
1452         simVertices.at(caloParticle.g4Tracks()[0].vertIndex()).position().z(), eta);
1453   }
1454 
1455   if (histograms.h_caloparticle_energy.count(pdgid)) {
1456     histograms.h_caloparticle_energy.at(pdgid)->Fill(caloParticle.energy());
1457   }
1458   if (histograms.h_caloparticle_pt.count(pdgid)) {
1459     histograms.h_caloparticle_pt.at(pdgid)->Fill(caloParticle.pt());
1460   }
1461   if (histograms.h_caloparticle_phi.count(pdgid)) {
1462     histograms.h_caloparticle_phi.at(pdgid)->Fill(caloParticle.phi());
1463   }
1464 
1465   if (histograms.h_caloparticle_nSimClusters.count(pdgid)) {
1466     histograms.h_caloparticle_nSimClusters.at(pdgid)->Fill(caloParticle.simClusters().size());
1467 
1468     int simHits = 0;
1469     int minLayerId = 999;
1470     int maxLayerId = 0;
1471 
1472     int simHits_matched = 0;
1473     int minLayerId_matched = 999;
1474     int maxLayerId_matched = 0;
1475 
1476     float energy = 0.;
1477     std::map<int, double> totenergy_layer;
1478 
1479     float hitEnergyWeight_invSum = 0;
1480     std::vector<std::pair<DetId, float>> haf_cp;
1481     for (const auto& sc : caloParticle.simClusters()) {
1482       LogDebug("HGCalValidator") << " This sim cluster has " << sc->hits_and_fractions().size() << " simHits and "
1483                                  << sc->energy() << " energy. " << std::endl;
1484       simHits += sc->hits_and_fractions().size();
1485       for (auto const& h_and_f : sc->hits_and_fractions()) {
1486         const auto hitDetId = h_and_f.first;
1487         const int layerId =
1488             recHitTools_->getLayerWithOffset(hitDetId) + layers * ((recHitTools_->zside(hitDetId) + 1) >> 1) - 1;
1489         // set to 0 if matched RecHit not found
1490         int layerId_matched_min = 999;
1491         int layerId_matched_max = 0;
1492         std::unordered_map<DetId, const HGCRecHit*>::const_iterator itcheck = hitMap.find(hitDetId);
1493         if (itcheck != hitMap.end()) {
1494           layerId_matched_min = layerId;
1495           layerId_matched_max = layerId;
1496           simHits_matched++;
1497 
1498           const auto hitEn = itcheck->second->energy();
1499           hitEnergyWeight_invSum += pow(hitEn, 2);
1500           const auto hitFr = h_and_f.second;
1501           const auto hitEnFr = hitEn * hitFr;
1502           energy += hitEnFr;
1503           histograms.h_caloparticle_nHits_matched_energy.at(pdgid)->Fill(hitEnFr);
1504           histograms.h_caloparticle_nHits_matched_energy_layer.at(pdgid)->Fill(layerId, hitEnFr);
1505 
1506           if (totenergy_layer.find(layerId) != totenergy_layer.end()) {
1507             totenergy_layer[layerId] = totenergy_layer.at(layerId) + hitEn;
1508           } else {
1509             totenergy_layer.emplace(layerId, hitEn);
1510           }
1511           if (caloParticle.simClusters().size() == 1)
1512             histograms.h_caloparticle_nHits_matched_energy_layer_1SimCl.at(pdgid)->Fill(layerId, hitEnFr);
1513 
1514           auto found = std::find_if(std::begin(haf_cp),
1515                                     std::end(haf_cp),
1516                                     [&hitDetId](const std::pair<DetId, float>& v) { return v.first == hitDetId; });
1517           if (found != haf_cp.end())
1518             found->second += hitFr;
1519           else
1520             haf_cp.emplace_back(hitDetId, hitFr);
1521 
1522         } else {
1523           LogDebug("HGCalValidator") << "   matched to RecHit NOT found !" << std::endl;
1524         }
1525 
1526         minLayerId = std::min(minLayerId, layerId);
1527         maxLayerId = std::max(maxLayerId, layerId);
1528         minLayerId_matched = std::min(minLayerId_matched, layerId_matched_min);
1529         maxLayerId_matched = std::max(maxLayerId_matched, layerId_matched_max);
1530       }
1531       LogDebug("HGCalValidator") << std::endl;
1532     }  // End loop over SimClusters of CaloParticle
1533     if (hitEnergyWeight_invSum)
1534       hitEnergyWeight_invSum = 1 / hitEnergyWeight_invSum;
1535 
1536     histograms.h_caloparticle_firstlayer.at(pdgid)->Fill(minLayerId);
1537     histograms.h_caloparticle_lastlayer.at(pdgid)->Fill(maxLayerId);
1538     histograms.h_caloparticle_layersnum.at(pdgid)->Fill(int(maxLayerId - minLayerId));
1539 
1540     histograms.h_caloparticle_firstlayer_matchedtoRecHit.at(pdgid)->Fill(minLayerId_matched);
1541     histograms.h_caloparticle_lastlayer_matchedtoRecHit.at(pdgid)->Fill(maxLayerId_matched);
1542     histograms.h_caloparticle_layersnum_matchedtoRecHit.at(pdgid)->Fill(int(maxLayerId_matched - minLayerId_matched));
1543 
1544     histograms.h_caloparticle_nHitsInSimClusters.at(pdgid)->Fill((float)simHits);
1545     histograms.h_caloparticle_nHitsInSimClusters_matchedtoRecHit.at(pdgid)->Fill((float)simHits_matched);
1546     histograms.h_caloparticle_selfenergy.at(pdgid)->Fill((float)energy);
1547     histograms.h_caloparticle_energyDifference.at(pdgid)->Fill((float)1. - energy / caloParticle.energy());
1548 
1549     //Calculate sum energy per-layer
1550     auto i = totenergy_layer.begin();
1551     double sum_energy = 0.0;
1552     while (i != totenergy_layer.end()) {
1553       sum_energy += i->second;
1554       histograms.h_caloparticle_sum_energy_layer.at(pdgid)->Fill(i->first, sum_energy / caloParticle.energy() * 100.);
1555       i++;
1556     }
1557 
1558     for (auto const& haf : haf_cp) {
1559       const auto hitEn = hitMap.find(haf.first)->second->energy();
1560       const auto weight = pow(hitEn, 2);
1561       histograms.h_caloparticle_fractions.at(pdgid)->Fill(haf.second, weight * hitEnergyWeight_invSum);
1562       histograms.h_caloparticle_fractions_weight.at(pdgid)->Fill(haf.second, weight * hitEnergyWeight_invSum, weight);
1563     }
1564   }
1565 }
1566 
1567 void HGVHistoProducerAlgo::HGVHistoProducerAlgo::fill_simCluster_histos(const Histograms& histograms,
1568                                                                         std::vector<SimCluster> const& simClusters,
1569                                                                         unsigned int layers,
1570                                                                         std::vector<int> thicknesses) const {
1571   //Each event to be treated as two events: an event in +ve endcap,
1572   //plus another event in -ve endcap. In this spirit there will be
1573   //a layer variable (layerid) that maps the layers in :
1574   //-z: 0->49
1575   //+z: 50->99
1576 
1577   //To keep track of total num of simClusters per layer
1578   //tnscpl[layerid]
1579   std::vector<int> tnscpl(1000, 0);  //tnscpl.clear(); tnscpl.reserve(1000);
1580 
1581   //To keep track of the total num of clusters per thickness in plus and in minus endcaps
1582   std::map<std::string, int> tnscpthplus;
1583   tnscpthplus.clear();
1584   std::map<std::string, int> tnscpthminus;
1585   tnscpthminus.clear();
1586   //At the beginning of the event all layers should be initialized to zero total clusters per thickness
1587   for (std::vector<int>::iterator it = thicknesses.begin(); it != thicknesses.end(); ++it) {
1588     tnscpthplus.insert(std::pair<std::string, int>(std::to_string(*it), 0));
1589     tnscpthminus.insert(std::pair<std::string, int>(std::to_string(*it), 0));
1590   }
1591   //To keep track of the total num of simClusters with mixed thickness hits per event
1592   tnscpthplus.insert(std::pair<std::string, int>("mixed", 0));
1593   tnscpthminus.insert(std::pair<std::string, int>("mixed", 0));
1594 
1595   //loop through simClusters
1596   for (const auto& sc : simClusters) {
1597     //Auxillary variables to count the number of different kind of hits in each simCluster
1598     int nthhits120p = 0;
1599     int nthhits200p = 0;
1600     int nthhits300p = 0;
1601     int nthhitsscintp = 0;
1602     int nthhits120m = 0;
1603     int nthhits200m = 0;
1604     int nthhits300m = 0;
1605     int nthhitsscintm = 0;
1606     //For the hits thickness of the layer cluster.
1607     double thickness = 0.;
1608     //To keep track if we added the simCluster in a specific layer
1609     std::vector<int> occurenceSCinlayer(1000, 0);  //[layerid][0 if not added]
1610 
1611     //loop through hits of the simCluster
1612     for (const auto& hAndF : sc.hits_and_fractions()) {
1613       const DetId sh_detid = hAndF.first;
1614 
1615       //The layer the cluster belongs to. As mentioned in the mapping above, it takes into account -z and +z.
1616       int layerid =
1617           recHitTools_->getLayerWithOffset(sh_detid) + layers * ((recHitTools_->zside(sh_detid) + 1) >> 1) - 1;
1618       //zside that the current cluster belongs to.
1619       int zside = recHitTools_->zside(sh_detid);
1620 
1621       //add the simCluster to the relevant layer. A SimCluster may give contribution to several layers.
1622       if (occurenceSCinlayer[layerid] == 0) {
1623         tnscpl[layerid]++;
1624       }
1625       occurenceSCinlayer[layerid]++;
1626 
1627       if (sh_detid.det() == DetId::Forward || sh_detid.det() == DetId::HGCalEE || sh_detid.det() == DetId::HGCalHSi) {
1628         thickness = recHitTools_->getSiThickness(sh_detid);
1629       } else if (sh_detid.det() == DetId::HGCalHSc) {
1630         thickness = -1;
1631       } else {
1632         LogDebug("HGCalValidator") << "These are HGCal simClusters, you shouldn't be here !!! " << layerid << "\n";
1633         continue;
1634       }
1635 
1636       if ((thickness == 120.) && (zside > 0.)) {
1637         nthhits120p++;
1638       } else if ((thickness == 120.) && (zside < 0.)) {
1639         nthhits120m++;
1640       } else if ((thickness == 200.) && (zside > 0.)) {
1641         nthhits200p++;
1642       } else if ((thickness == 200.) && (zside < 0.)) {
1643         nthhits200m++;
1644       } else if ((thickness == 300.) && (zside > 0.)) {
1645         nthhits300p++;
1646       } else if ((thickness == 300.) && (zside < 0.)) {
1647         nthhits300m++;
1648       } else if ((thickness == -1) && (zside > 0.)) {
1649         nthhitsscintp++;
1650       } else if ((thickness == -1) && (zside < 0.)) {
1651         nthhitsscintm++;
1652       } else {  //assert(0);
1653         LogDebug("HGCalValidator")
1654             << " You are running a geometry that contains thicknesses different than the normal ones. "
1655             << "\n";
1656       }
1657 
1658     }  //end of loop through hits
1659 
1660     //Check for simultaneously having hits of different kind. Checking at least two combinations is sufficient.
1661     if ((nthhits120p != 0 && nthhits200p != 0) || (nthhits120p != 0 && nthhits300p != 0) ||
1662         (nthhits120p != 0 && nthhitsscintp != 0) || (nthhits200p != 0 && nthhits300p != 0) ||
1663         (nthhits200p != 0 && nthhitsscintp != 0) || (nthhits300p != 0 && nthhitsscintp != 0)) {
1664       tnscpthplus["mixed"]++;
1665     } else if ((nthhits120p != 0 || nthhits200p != 0 || nthhits300p != 0 || nthhitsscintp != 0)) {
1666       //This is a cluster with hits of one kind
1667       tnscpthplus[std::to_string((int)thickness)]++;
1668     }
1669     if ((nthhits120m != 0 && nthhits200m != 0) || (nthhits120m != 0 && nthhits300m != 0) ||
1670         (nthhits120m != 0 && nthhitsscintm != 0) || (nthhits200m != 0 && nthhits300m != 0) ||
1671         (nthhits200m != 0 && nthhitsscintm != 0) || (nthhits300m != 0 && nthhitsscintm != 0)) {
1672       tnscpthminus["mixed"]++;
1673     } else if ((nthhits120m != 0 || nthhits200m != 0 || nthhits300m != 0 || nthhitsscintm != 0)) {
1674       //This is a cluster with hits of one kind
1675       tnscpthminus[std::to_string((int)thickness)]++;
1676     }
1677 
1678   }  //end of loop through SimClusters of the event
1679 
1680   //Per layer : Loop 0->99
1681   for (unsigned ilayer = 0; ilayer < layers * 2; ++ilayer)
1682     if (histograms.h_simclusternum_perlayer.count(ilayer))
1683       histograms.h_simclusternum_perlayer.at(ilayer)->Fill(tnscpl[ilayer]);
1684 
1685   //Per thickness
1686   for (std::vector<int>::iterator it = thicknesses.begin(); it != thicknesses.end(); ++it) {
1687     if (histograms.h_simclusternum_perthick.count(*it)) {
1688       histograms.h_simclusternum_perthick.at(*it)->Fill(tnscpthplus[std::to_string(*it)]);
1689       histograms.h_simclusternum_perthick.at(*it)->Fill(tnscpthminus[std::to_string(*it)]);
1690     }
1691   }
1692   //Mixed thickness clusters
1693   histograms.h_mixedhitssimcluster_zplus->Fill(tnscpthplus["mixed"]);
1694   histograms.h_mixedhitssimcluster_zminus->Fill(tnscpthminus["mixed"]);
1695 }
1696 
1697 void HGVHistoProducerAlgo::HGVHistoProducerAlgo::fill_simClusterAssociation_histos(
1698     const Histograms& histograms,
1699     const int count,
1700     edm::Handle<reco::CaloClusterCollection> clusterHandle,
1701     const reco::CaloClusterCollection& clusters,
1702     edm::Handle<std::vector<SimCluster>> simClusterHandle,
1703     std::vector<SimCluster> const& simClusters,
1704     std::vector<size_t> const& sCIndices,
1705     const std::vector<float>& mask,
1706     std::unordered_map<DetId, const HGCRecHit*> const& hitMap,
1707     unsigned int layers,
1708     const hgcal::RecoToSimCollectionWithSimClusters& scsInLayerClusterMap,
1709     const hgcal::SimToRecoCollectionWithSimClusters& lcsInSimClusterMap) const {
1710   //Each event to be treated as two events: an event in +ve endcap,
1711   //plus another event in -ve endcap. In this spirit there will be
1712   //a layer variable (layerid) that maps the layers in :
1713   //-z: 0->49
1714   //+z: 50->99
1715 
1716   //Will add some general plots on the specific mask in the future.
1717 
1718   layerClusters_to_SimClusters(histograms,
1719                                count,
1720                                clusterHandle,
1721                                clusters,
1722                                simClusterHandle,
1723                                simClusters,
1724                                sCIndices,
1725                                mask,
1726                                hitMap,
1727                                layers,
1728                                scsInLayerClusterMap,
1729                                lcsInSimClusterMap);
1730 }
1731 
1732 void HGVHistoProducerAlgo::fill_cluster_histos(const Histograms& histograms,
1733                                                const int count,
1734                                                const reco::CaloCluster& cluster) const {
1735   const auto eta = getEta(cluster.eta());
1736   histograms.h_cluster_eta[count]->Fill(eta);
1737 }
1738 
1739 void HGVHistoProducerAlgo::layerClusters_to_CaloParticles(const Histograms& histograms,
1740                                                           edm::Handle<reco::CaloClusterCollection> clusterHandle,
1741                                                           const reco::CaloClusterCollection& clusters,
1742                                                           edm::Handle<std::vector<CaloParticle>> caloParticleHandle,
1743                                                           std::vector<CaloParticle> const& cP,
1744                                                           std::vector<size_t> const& cPIndices,
1745                                                           std::vector<size_t> const& cPSelectedIndices,
1746                                                           std::unordered_map<DetId, const HGCRecHit*> const& hitMap,
1747                                                           unsigned int layers,
1748                                                           const hgcal::RecoToSimCollection& cpsInLayerClusterMap,
1749                                                           const hgcal::SimToRecoCollection& cPOnLayerMap) const {
1750   const auto nLayerClusters = clusters.size();
1751 
1752   std::unordered_map<DetId, std::vector<HGVHistoProducerAlgo::detIdInfoInCluster>> detIdToCaloParticleId_Map;
1753   std::unordered_map<DetId, std::vector<HGVHistoProducerAlgo::detIdInfoInCluster>> detIdToLayerClusterId_Map;
1754 
1755   // The association has to be done in an all-vs-all fashion.
1756   // For this reason use the full set of CaloParticles, with the only filter on bx
1757   for (const auto& cpId : cPIndices) {
1758     for (const auto& simCluster : cP[cpId].simClusters()) {
1759       for (const auto& it_haf : simCluster->hits_and_fractions()) {
1760         const DetId hitid = (it_haf.first);
1761         if (hitMap.find(hitid) != hitMap.end()) {
1762           if (detIdToCaloParticleId_Map.find(hitid) == detIdToCaloParticleId_Map.end()) {
1763             detIdToCaloParticleId_Map[hitid] = std::vector<HGVHistoProducerAlgo::detIdInfoInCluster>();
1764             detIdToCaloParticleId_Map[hitid].emplace_back(
1765                 HGVHistoProducerAlgo::detIdInfoInCluster{cpId, it_haf.second});
1766           } else {
1767             auto findHitIt =
1768                 std::find(detIdToCaloParticleId_Map[hitid].begin(),
1769                           detIdToCaloParticleId_Map[hitid].end(),
1770                           HGVHistoProducerAlgo::detIdInfoInCluster{
1771                               cpId, 0.f});  // only the first element is used for the matching (overloaded operator==)
1772             if (findHitIt != detIdToCaloParticleId_Map[hitid].end())
1773               findHitIt->fraction += it_haf.second;
1774             else
1775               detIdToCaloParticleId_Map[hitid].emplace_back(
1776                   HGVHistoProducerAlgo::detIdInfoInCluster{cpId, it_haf.second});
1777           }
1778         }
1779       }
1780     }
1781   }
1782 
1783   for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) {
1784     const auto& hits_and_fractions = clusters[lcId].hitsAndFractions();
1785     const auto numberOfHitsInLC = hits_and_fractions.size();
1786 
1787     // This vector will store, for each hit in the Layercluster, the index of
1788     // the CaloParticle that contributed the most, in terms of energy, to it.
1789     // Special values are:
1790     //
1791     // -2  --> the reconstruction fraction of the RecHit is 0 (used in the past to monitor Halo Hits)
1792     // -3  --> same as before with the added condition that no CaloParticle has been linked to this RecHit
1793     // -1  --> the reco fraction is >0, but no CaloParticle has been linked to it
1794     // >=0 --> index of the linked CaloParticle
1795     std::vector<int> hitsToCaloParticleId(numberOfHitsInLC);
1796     const auto firstHitDetId = hits_and_fractions[0].first;
1797     int lcLayerId =
1798         recHitTools_->getLayerWithOffset(firstHitDetId) + layers * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
1799 
1800     // This will store the fraction of the CaloParticle energy shared with the LayerCluster: e_shared/cp_energy
1801     std::unordered_map<unsigned, float> CPEnergyInLC;
1802 
1803     for (unsigned int iHit = 0; iHit < numberOfHitsInLC; iHit++) {
1804       const DetId rh_detid = hits_and_fractions[iHit].first;
1805       const auto rhFraction = hits_and_fractions[iHit].second;
1806 
1807       std::unordered_map<DetId, const HGCRecHit*>::const_iterator itcheck = hitMap.find(rh_detid);
1808       const HGCRecHit* hit = itcheck->second;
1809 
1810       if (detIdToLayerClusterId_Map.find(rh_detid) == detIdToLayerClusterId_Map.end()) {
1811         detIdToLayerClusterId_Map[rh_detid] = std::vector<HGVHistoProducerAlgo::detIdInfoInCluster>();
1812       }
1813       detIdToLayerClusterId_Map[rh_detid].emplace_back(HGVHistoProducerAlgo::detIdInfoInCluster{lcId, rhFraction});
1814 
1815       const auto& hit_find_in_CP = detIdToCaloParticleId_Map.find(rh_detid);
1816 
1817       // if the fraction is zero or the hit does not belong to any calo
1818       // particle, set the caloparticleId for the hit to -1 this will
1819       // contribute to the number of noise hits
1820 
1821       // MR Remove the case in which the fraction is 0, since this could be a
1822       // real hit that has been marked as halo.
1823       if (rhFraction == 0.) {
1824         hitsToCaloParticleId[iHit] = -2;
1825       }
1826       if (hit_find_in_CP == detIdToCaloParticleId_Map.end()) {
1827         hitsToCaloParticleId[iHit] -= 1;
1828       } else {
1829         auto maxCPEnergyInLC = 0.f;
1830         auto maxCPId = -1;
1831         for (auto& h : hit_find_in_CP->second) {
1832           const auto iCP = h.clusterId;
1833           CPEnergyInLC[iCP] += h.fraction * hit->energy();
1834           // Keep track of which CaloParticle contributed the most, in terms
1835           // of energy, to this specific LayerCluster.
1836           if (CPEnergyInLC[iCP] > maxCPEnergyInLC) {
1837             maxCPEnergyInLC = CPEnergyInLC[iCP];
1838             maxCPId = iCP;
1839           }
1840         }
1841         hitsToCaloParticleId[iHit] = maxCPId;
1842       }
1843       histograms.h_cellAssociation_perlayer.at(lcLayerId)->Fill(
1844           hitsToCaloParticleId[iHit] > 0. ? 0. : hitsToCaloParticleId[iHit]);
1845     }  // End loop over hits on a LayerCluster
1846 
1847   }  // End of loop over LayerClusters
1848 
1849   // Fill the plots to compute the different metrics linked to
1850   // reco-level, namely fake-rate an merge-rate. In this loop should *not*
1851   // restrict only to the selected caloParaticles.
1852   for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) {
1853     const auto firstHitDetId = (clusters[lcId].hitsAndFractions())[0].first;
1854     const int lcLayerId =
1855         recHitTools_->getLayerWithOffset(firstHitDetId) + layers * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
1856     histograms.h_denom_layercl_eta_perlayer.at(lcLayerId)->Fill(clusters[lcId].eta());
1857     histograms.h_denom_layercl_phi_perlayer.at(lcLayerId)->Fill(clusters[lcId].phi());
1858     //
1859     const edm::Ref<reco::CaloClusterCollection> lcRef(clusterHandle, lcId);
1860     const auto& cpsIt = cpsInLayerClusterMap.find(lcRef);
1861     if (cpsIt == cpsInLayerClusterMap.end())
1862       continue;
1863 
1864     const auto lc_en = clusters[lcId].energy();
1865     const auto& cps = cpsIt->val;
1866     if (lc_en == 0. && !cps.empty()) {
1867       for (const auto& cpPair : cps)
1868         histograms.h_score_layercl2caloparticle_perlayer.at(lcLayerId)->Fill(cpPair.second);
1869       continue;
1870     }
1871     for (const auto& cpPair : cps) {
1872       LogDebug("HGCalValidator") << "layerCluster Id: \t" << lcId << "\t CP id: \t" << cpPair.first.index()
1873                                  << "\t score \t" << cpPair.second << std::endl;
1874       histograms.h_score_layercl2caloparticle_perlayer.at(lcLayerId)->Fill(cpPair.second);
1875       auto const& cp_linked =
1876           std::find_if(std::begin(cPOnLayerMap[cpPair.first]),
1877                        std::end(cPOnLayerMap[cpPair.first]),
1878                        [&lcRef](const std::pair<edm::Ref<reco::CaloClusterCollection>, std::pair<float, float>>& p) {
1879                          return p.first == lcRef;
1880                        });
1881       if (cp_linked ==
1882           cPOnLayerMap[cpPair.first].end())  // This should never happen by construction of the association maps
1883         continue;
1884       histograms.h_sharedenergy_layercl2caloparticle_perlayer.at(lcLayerId)->Fill(cp_linked->second.first / lc_en,
1885                                                                                   lc_en);
1886       histograms.h_energy_vs_score_layercl2caloparticle_perlayer.at(lcLayerId)->Fill(cpPair.second,
1887                                                                                      cp_linked->second.first / lc_en);
1888     }
1889     const auto assoc =
1890         std::count_if(std::begin(cps), std::end(cps), [](const auto& obj) { return obj.second < ScoreCutLCtoCP_; });
1891     if (assoc) {
1892       histograms.h_num_layercl_eta_perlayer.at(lcLayerId)->Fill(clusters[lcId].eta());
1893       histograms.h_num_layercl_phi_perlayer.at(lcLayerId)->Fill(clusters[lcId].phi());
1894       if (assoc > 1) {
1895         histograms.h_numMerge_layercl_eta_perlayer.at(lcLayerId)->Fill(clusters[lcId].eta());
1896         histograms.h_numMerge_layercl_phi_perlayer.at(lcLayerId)->Fill(clusters[lcId].phi());
1897       }
1898       const auto& best = std::min_element(
1899           std::begin(cps), std::end(cps), [](const auto& obj1, const auto& obj2) { return obj1.second < obj2.second; });
1900       const auto& best_cp_linked =
1901           std::find_if(std::begin(cPOnLayerMap[best->first]),
1902                        std::end(cPOnLayerMap[best->first]),
1903                        [&lcRef](const std::pair<edm::Ref<reco::CaloClusterCollection>, std::pair<float, float>>& p) {
1904                          return p.first == lcRef;
1905                        });
1906       if (best_cp_linked ==
1907           cPOnLayerMap[best->first].end())  // This should never happen by construction of the association maps
1908         continue;
1909       histograms.h_sharedenergy_layercl2caloparticle_vs_eta_perlayer.at(lcLayerId)->Fill(
1910           clusters[lcId].eta(), best_cp_linked->second.first / lc_en);
1911       histograms.h_sharedenergy_layercl2caloparticle_vs_phi_perlayer.at(lcLayerId)->Fill(
1912           clusters[lcId].phi(), best_cp_linked->second.first / lc_en);
1913     }
1914   }  // End of loop over LayerClusters
1915 
1916   // Here Fill the plots to compute the different metrics linked to
1917   // gen-level, namely efficiency and duplicate. In this loop should restrict
1918   // only to the selected caloParaticles.
1919   for (const auto& cpId : cPSelectedIndices) {
1920     const edm::Ref<CaloParticleCollection> cpRef(caloParticleHandle, cpId);
1921     const auto& lcsIt = cPOnLayerMap.find(cpRef);
1922 
1923     std::map<unsigned int, float> cPEnergyOnLayer;
1924     for (unsigned int layerId = 0; layerId < layers * 2; ++layerId)
1925       cPEnergyOnLayer[layerId] = 0;
1926 
1927     for (const auto& simCluster : cP[cpId].simClusters()) {
1928       for (const auto& it_haf : simCluster->hits_and_fractions()) {
1929         const DetId hitid = (it_haf.first);
1930         const auto hitLayerId =
1931             recHitTools_->getLayerWithOffset(hitid) + layers * ((recHitTools_->zside(hitid) + 1) >> 1) - 1;
1932         std::unordered_map<DetId, const HGCRecHit*>::const_iterator itcheck = hitMap.find(hitid);
1933         if (itcheck != hitMap.end()) {
1934           const HGCRecHit* hit = itcheck->second;
1935           cPEnergyOnLayer[hitLayerId] += it_haf.second * hit->energy();
1936         }
1937       }
1938     }
1939 
1940     for (unsigned int layerId = 0; layerId < layers * 2; ++layerId) {
1941       if (!cPEnergyOnLayer[layerId])
1942         continue;
1943 
1944       histograms.h_denom_caloparticle_eta_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().eta());
1945       histograms.h_denom_caloparticle_phi_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().phi());
1946 
1947       if (lcsIt == cPOnLayerMap.end())
1948         continue;
1949       const auto& lcs = lcsIt->val;
1950 
1951       auto getLCLayerId = [&](const unsigned int lcId) {
1952         const auto firstHitDetId = (clusters[lcId].hitsAndFractions())[0].first;
1953         const auto lcLayerId = recHitTools_->getLayerWithOffset(firstHitDetId) +
1954                                layers * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
1955         return lcLayerId;
1956       };
1957 
1958       for (const auto& lcPair : lcs) {
1959         if (getLCLayerId(lcPair.first.index()) != layerId)
1960           continue;
1961         histograms.h_score_caloparticle2layercl_perlayer.at(layerId)->Fill(lcPair.second.second);
1962         histograms.h_sharedenergy_caloparticle2layercl_perlayer.at(layerId)->Fill(
1963             lcPair.second.first / cPEnergyOnLayer[layerId], cPEnergyOnLayer[layerId]);
1964         histograms.h_energy_vs_score_caloparticle2layercl_perlayer.at(layerId)->Fill(
1965             lcPair.second.second, lcPair.second.first / cPEnergyOnLayer[layerId]);
1966       }
1967       const auto assoc = std::count_if(std::begin(lcs), std::end(lcs), [&](const auto& obj) {
1968         if (getLCLayerId(obj.first.index()) != layerId)
1969           return false;
1970         else
1971           return obj.second.second < ScoreCutCPtoLC_;
1972       });
1973       if (assoc) {
1974         histograms.h_num_caloparticle_eta_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().eta());
1975         histograms.h_num_caloparticle_phi_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().phi());
1976         if (assoc > 1) {
1977           histograms.h_numDup_caloparticle_eta_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().eta());
1978           histograms.h_numDup_caloparticle_phi_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().phi());
1979         }
1980         const auto best = std::min_element(std::begin(lcs), std::end(lcs), [&](const auto& obj1, const auto& obj2) {
1981           if (getLCLayerId(obj1.first.index()) != layerId)
1982             return false;
1983           else if (getLCLayerId(obj2.first.index()) == layerId)
1984             return obj1.second.second < obj2.second.second;
1985           else
1986             return true;
1987         });
1988         histograms.h_sharedenergy_caloparticle2layercl_vs_eta_perlayer.at(layerId)->Fill(
1989             cP[cpId].g4Tracks()[0].momentum().eta(), best->second.first / cPEnergyOnLayer[layerId]);
1990         histograms.h_sharedenergy_caloparticle2layercl_vs_phi_perlayer.at(layerId)->Fill(
1991             cP[cpId].g4Tracks()[0].momentum().phi(), best->second.first / cPEnergyOnLayer[layerId]);
1992       }
1993     }
1994   }
1995 }
1996 
1997 void HGVHistoProducerAlgo::layerClusters_to_SimClusters(
1998     const Histograms& histograms,
1999     const int count,
2000     edm::Handle<reco::CaloClusterCollection> clusterHandle,
2001     const reco::CaloClusterCollection& clusters,
2002     edm::Handle<std::vector<SimCluster>> simClusterHandle,
2003     std::vector<SimCluster> const& sC,
2004     std::vector<size_t> const& sCIndices,
2005     const std::vector<float>& mask,
2006     std::unordered_map<DetId, const HGCRecHit*> const& hitMap,
2007     unsigned int layers,
2008     const hgcal::RecoToSimCollectionWithSimClusters& scsInLayerClusterMap,
2009     const hgcal::SimToRecoCollectionWithSimClusters& lcsInSimClusterMap) const {
2010   // Here fill the plots to compute the different metrics linked to
2011   // reco-level, namely fake-rate and merge-rate. In this loop should *not*
2012   // restrict only to the selected SimClusters.
2013   for (unsigned int lcId = 0; lcId < clusters.size(); ++lcId) {
2014     if (mask[lcId] != 0.) {
2015       LogDebug("HGCalValidator") << "Skipping layer cluster " << lcId << " not belonging to mask" << std::endl;
2016       continue;
2017     }
2018     const auto firstHitDetId = (clusters[lcId].hitsAndFractions())[0].first;
2019     const auto lcLayerId =
2020         recHitTools_->getLayerWithOffset(firstHitDetId) + layers * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
2021     //Although the ones below are already created in the LC to CP association, will
2022     //recreate them here since in the post processor it looks in a specific directory.
2023     histograms.h_denom_layercl_in_simcl_eta_perlayer[count].at(lcLayerId)->Fill(clusters[lcId].eta());
2024     histograms.h_denom_layercl_in_simcl_phi_perlayer[count].at(lcLayerId)->Fill(clusters[lcId].phi());
2025     //
2026     const edm::Ref<reco::CaloClusterCollection> lcRef(clusterHandle, lcId);
2027     const auto& scsIt = scsInLayerClusterMap.find(lcRef);
2028     if (scsIt == scsInLayerClusterMap.end())
2029       continue;
2030 
2031     const auto lc_en = clusters[lcId].energy();
2032     const auto& scs = scsIt->val;
2033     // If a reconstructed LayerCluster has energy 0 but is linked to at least a
2034     // SimCluster, then his score should be 1 as set in the associator
2035     if (lc_en == 0. && !scs.empty()) {
2036       for (const auto& scPair : scs) {
2037         histograms.h_score_layercl2simcluster_perlayer[count].at(lcLayerId)->Fill(scPair.second);
2038       }
2039       continue;
2040     }
2041     //Loop through all SimClusters linked to the layer cluster under study
2042     for (const auto& scPair : scs) {
2043       LogDebug("HGCalValidator") << "layerCluster Id: \t" << lcId << "\t SC id: \t" << scPair.first.index()
2044                                  << "\t score \t" << scPair.second << std::endl;
2045       //This should be filled #layerClusters in layer x #linked SimClusters
2046       histograms.h_score_layercl2simcluster_perlayer[count].at(lcLayerId)->Fill(scPair.second);
2047       auto const& sc_linked =
2048           std::find_if(std::begin(lcsInSimClusterMap[scPair.first]),
2049                        std::end(lcsInSimClusterMap[scPair.first]),
2050                        [&lcRef](const std::pair<edm::Ref<reco::CaloClusterCollection>, std::pair<float, float>>& p) {
2051                          return p.first == lcRef;
2052                        });
2053       if (sc_linked ==
2054           lcsInSimClusterMap[scPair.first].end())  // This should never happen by construction of the association maps
2055         continue;
2056       histograms.h_sharedenergy_layercl2simcluster_perlayer[count].at(lcLayerId)->Fill(sc_linked->second.first / lc_en,
2057                                                                                        lc_en);
2058       histograms.h_energy_vs_score_layercl2simcluster_perlayer[count].at(lcLayerId)->Fill(
2059           scPair.second, sc_linked->second.first / lc_en);
2060     }
2061     //Here he counts how many of the linked SimClusters of the layer cluster under study have a score above a certain value.
2062     const auto assoc =
2063         std::count_if(std::begin(scs), std::end(scs), [](const auto& obj) { return obj.second < ScoreCutLCtoSC_; });
2064     if (assoc) {
2065       histograms.h_num_layercl_in_simcl_eta_perlayer[count].at(lcLayerId)->Fill(clusters[lcId].eta());
2066       histograms.h_num_layercl_in_simcl_phi_perlayer[count].at(lcLayerId)->Fill(clusters[lcId].phi());
2067       if (assoc > 1) {
2068         histograms.h_numMerge_layercl_in_simcl_eta_perlayer[count].at(lcLayerId)->Fill(clusters[lcId].eta());
2069         histograms.h_numMerge_layercl_in_simcl_phi_perlayer[count].at(lcLayerId)->Fill(clusters[lcId].phi());
2070       }
2071       const auto& best = std::min_element(
2072           std::begin(scs), std::end(scs), [](const auto& obj1, const auto& obj2) { return obj1.second < obj2.second; });
2073       //From all SimClusters he founds the one with the best (lowest) score and takes his scId
2074       const auto& best_sc_linked =
2075           std::find_if(std::begin(lcsInSimClusterMap[best->first]),
2076                        std::end(lcsInSimClusterMap[best->first]),
2077                        [&lcRef](const std::pair<edm::Ref<reco::CaloClusterCollection>, std::pair<float, float>>& p) {
2078                          return p.first == lcRef;
2079                        });
2080       if (best_sc_linked ==
2081           lcsInSimClusterMap[best->first].end())  // This should never happen by construction of the association maps
2082         continue;
2083       histograms.h_sharedenergy_layercl2simcluster_vs_eta_perlayer[count].at(lcLayerId)->Fill(
2084           clusters[lcId].eta(), best_sc_linked->second.first / lc_en);
2085       histograms.h_sharedenergy_layercl2simcluster_vs_phi_perlayer[count].at(lcLayerId)->Fill(
2086           clusters[lcId].phi(), best_sc_linked->second.first / lc_en);
2087     }
2088   }  // End of loop over LayerClusters
2089 
2090   // Fill the plots to compute the different metrics linked to
2091   // gen-level, namely efficiency and duplicate. In this loop should restrict
2092   // only to the selected SimClusters.
2093   for (const auto& scId : sCIndices) {
2094     const edm::Ref<SimClusterCollection> scRef(simClusterHandle, scId);
2095     const auto& lcsIt = lcsInSimClusterMap.find(scRef);
2096 
2097     std::map<unsigned int, float> sCEnergyOnLayer;
2098     for (unsigned int layerId = 0; layerId < layers * 2; ++layerId)
2099       sCEnergyOnLayer[layerId] = 0;
2100 
2101     for (const auto& it_haf : sC[scId].hits_and_fractions()) {
2102       const DetId hitid = (it_haf.first);
2103       const auto scLayerId =
2104           recHitTools_->getLayerWithOffset(hitid) + layers * ((recHitTools_->zside(hitid) + 1) >> 1) - 1;
2105       std::unordered_map<DetId, const HGCRecHit*>::const_iterator itcheck = hitMap.find(hitid);
2106       if (itcheck != hitMap.end()) {
2107         const HGCRecHit* hit = itcheck->second;
2108         sCEnergyOnLayer[scLayerId] += it_haf.second * hit->energy();
2109       }
2110     }
2111 
2112     for (unsigned int layerId = 0; layerId < layers * 2; ++layerId) {
2113       if (!sCEnergyOnLayer[layerId])
2114         continue;
2115 
2116       histograms.h_denom_simcluster_eta_perlayer[count].at(layerId)->Fill(sC[scId].eta());
2117       histograms.h_denom_simcluster_phi_perlayer[count].at(layerId)->Fill(sC[scId].phi());
2118 
2119       if (lcsIt == lcsInSimClusterMap.end())
2120         continue;
2121       const auto& lcs = lcsIt->val;
2122 
2123       auto getLCLayerId = [&](const unsigned int lcId) {
2124         const auto firstHitDetId = (clusters[lcId].hitsAndFractions())[0].first;
2125         const unsigned int lcLayerId = recHitTools_->getLayerWithOffset(firstHitDetId) +
2126                                        layers * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
2127         return lcLayerId;
2128       };
2129 
2130       //Loop through layer clusters linked to the SimCluster under study
2131       for (const auto& lcPair : lcs) {
2132         auto lcId = lcPair.first.index();
2133         if (mask[lcId] != 0.) {
2134           LogDebug("HGCalValidator") << "Skipping layer cluster " << lcId << " not belonging to mask" << std::endl;
2135           continue;
2136         }
2137 
2138         if (getLCLayerId(lcId) != layerId)
2139           continue;
2140         histograms.h_score_simcluster2layercl_perlayer[count].at(layerId)->Fill(lcPair.second.second);
2141         histograms.h_sharedenergy_simcluster2layercl_perlayer[count].at(layerId)->Fill(
2142             lcPair.second.first / sCEnergyOnLayer[layerId], sCEnergyOnLayer[layerId]);
2143         histograms.h_energy_vs_score_simcluster2layercl_perlayer[count].at(layerId)->Fill(
2144             lcPair.second.second, lcPair.second.first / sCEnergyOnLayer[layerId]);
2145       }
2146       const auto assoc = std::count_if(std::begin(lcs), std::end(lcs), [&](const auto& obj) {
2147         if (getLCLayerId(obj.first.index()) != layerId)
2148           return false;
2149         else
2150           return obj.second.second < ScoreCutSCtoLC_;
2151       });
2152       if (assoc) {
2153         histograms.h_num_simcluster_eta_perlayer[count].at(layerId)->Fill(sC[scId].eta());
2154         histograms.h_num_simcluster_phi_perlayer[count].at(layerId)->Fill(sC[scId].phi());
2155         if (assoc > 1) {
2156           histograms.h_numDup_simcluster_eta_perlayer[count].at(layerId)->Fill(sC[scId].eta());
2157           histograms.h_numDup_simcluster_phi_perlayer[count].at(layerId)->Fill(sC[scId].phi());
2158         }
2159         const auto best = std::min_element(std::begin(lcs), std::end(lcs), [&](const auto& obj1, const auto& obj2) {
2160           if (getLCLayerId(obj1.first.index()) != layerId)
2161             return false;
2162           else if (getLCLayerId(obj2.first.index()) == layerId)
2163             return obj1.second.second < obj2.second.second;
2164           else
2165             return true;
2166         });
2167         histograms.h_sharedenergy_simcluster2layercl_vs_eta_perlayer[count].at(layerId)->Fill(
2168             sC[scId].eta(), best->second.first / sCEnergyOnLayer[layerId]);
2169         histograms.h_sharedenergy_simcluster2layercl_vs_phi_perlayer[count].at(layerId)->Fill(
2170             sC[scId].phi(), best->second.first / sCEnergyOnLayer[layerId]);
2171       }
2172     }
2173   }
2174 }
2175 
2176 void HGVHistoProducerAlgo::fill_generic_cluster_histos(const Histograms& histograms,
2177                                                        const int count,
2178                                                        edm::Handle<reco::CaloClusterCollection> clusterHandle,
2179                                                        const reco::CaloClusterCollection& clusters,
2180                                                        const Density& densities,
2181                                                        edm::Handle<std::vector<CaloParticle>> caloParticleHandle,
2182                                                        std::vector<CaloParticle> const& cP,
2183                                                        std::vector<size_t> const& cPIndices,
2184                                                        std::vector<size_t> const& cPSelectedIndices,
2185                                                        std::unordered_map<DetId, const HGCRecHit*> const& hitMap,
2186                                                        std::map<double, double> cummatbudg,
2187                                                        unsigned int layers,
2188                                                        std::vector<int> thicknesses,
2189                                                        const hgcal::RecoToSimCollection& cpsInLayerClusterMap,
2190                                                        const hgcal::SimToRecoCollection& cPOnLayerMap) const {
2191   //Each event to be treated as two events: an event in +ve endcap,
2192   //plus another event in -ve endcap. In this spirit there will be
2193   //a layer variable (layerid) that maps the layers in :
2194   //-z: 0->51
2195   //+z: 52->103
2196 
2197   //To keep track of total num of layer clusters per layer
2198   //tnlcpl[layerid]
2199   std::vector<int> tnlcpl(1000, 0);  //tnlcpl.clear(); tnlcpl.reserve(1000);
2200 
2201   //To keep track of the total num of clusters per thickness in plus and in minus endcaps
2202   std::map<std::string, int> tnlcpthplus;
2203   tnlcpthplus.clear();
2204   std::map<std::string, int> tnlcpthminus;
2205   tnlcpthminus.clear();
2206   //At the beginning of the event all layers should be initialized to zero total clusters per thickness
2207   for (std::vector<int>::iterator it = thicknesses.begin(); it != thicknesses.end(); ++it) {
2208     tnlcpthplus.insert(std::pair<std::string, int>(std::to_string(*it), 0));
2209     tnlcpthminus.insert(std::pair<std::string, int>(std::to_string(*it), 0));
2210   }
2211   //To keep track of the total num of clusters with mixed thickness hits per event
2212   tnlcpthplus.insert(std::pair<std::string, int>("mixed", 0));
2213   tnlcpthminus.insert(std::pair<std::string, int>("mixed", 0));
2214 
2215   layerClusters_to_CaloParticles(histograms,
2216                                  clusterHandle,
2217                                  clusters,
2218                                  caloParticleHandle,
2219                                  cP,
2220                                  cPIndices,
2221                                  cPSelectedIndices,
2222                                  hitMap,
2223                                  layers,
2224                                  cpsInLayerClusterMap,
2225                                  cPOnLayerMap);
2226 
2227   //To find out the total amount of energy clustered per layer
2228   //Initialize with zeros because I see clear gives weird numbers.
2229   std::vector<double> tecpl(1000, 0.0);  //tecpl.clear(); tecpl.reserve(1000);
2230   //for the longitudinal depth barycenter
2231   std::vector<double> ldbar(1000, 0.0);  //ldbar.clear(); ldbar.reserve(1000);
2232 
2233   // Need to compare with the total amount of energy coming from CaloParticles
2234   double caloparteneplus = 0.;
2235   double caloparteneminus = 0.;
2236   for (const auto& cpId : cPIndices) {
2237     if (cP[cpId].eta() >= 0.) {
2238       caloparteneplus = caloparteneplus + cP[cpId].energy();
2239     } else if (cP[cpId].eta() < 0.) {
2240       caloparteneminus = caloparteneminus + cP[cpId].energy();
2241     }
2242   }
2243 
2244   // loop through clusters of the event
2245   for (const auto& lcId : clusters) {
2246     const auto seedid = lcId.seed();
2247     const double seedx = recHitTools_->getPosition(seedid).x();
2248     const double seedy = recHitTools_->getPosition(seedid).y();
2249     DetId maxid = findmaxhit(lcId, hitMap);
2250 
2251     // const DetId maxid = lcId.max();
2252     double maxx = recHitTools_->getPosition(maxid).x();
2253     double maxy = recHitTools_->getPosition(maxid).y();
2254 
2255     //Auxillary variables to count the number of different kind of hits in each cluster
2256     int nthhits120p = 0;
2257     int nthhits200p = 0;
2258     int nthhits300p = 0;
2259     int nthhitsscintp = 0;
2260     int nthhits120m = 0;
2261     int nthhits200m = 0;
2262     int nthhits300m = 0;
2263     int nthhitsscintm = 0;
2264     //For the hits thickness of the layer cluster.
2265     double thickness = 0.;
2266     //The layer the cluster belongs to. As mentioned in the mapping above, it takes into account -z and +z.
2267     int layerid = 0;
2268     // Need another layer variable for the longitudinal material budget file reading.
2269     //In this case need no distinction between -z and +z.
2270     int lay = 0;
2271     // Need to save the combination thick_lay
2272     std::string istr = "";
2273     //boolean to check for the layer that the cluster belong to. Maybe later will check all the layer hits.
2274     bool cluslay = true;
2275     //zside that the current cluster belongs to.
2276     int zside = 0;
2277 
2278     const auto& hits_and_fractions = lcId.hitsAndFractions();
2279     for (std::vector<std::pair<DetId, float>>::const_iterator it_haf = hits_and_fractions.begin();
2280          it_haf != hits_and_fractions.end();
2281          ++it_haf) {
2282       const DetId rh_detid = it_haf->first;
2283       //The layer that the current hit belongs to
2284       layerid = recHitTools_->getLayerWithOffset(rh_detid) + layers * ((recHitTools_->zside(rh_detid) + 1) >> 1) - 1;
2285       lay = recHitTools_->getLayerWithOffset(rh_detid);
2286       zside = recHitTools_->zside(rh_detid);
2287       if (rh_detid.det() == DetId::Forward || rh_detid.det() == DetId::HGCalEE || rh_detid.det() == DetId::HGCalHSi) {
2288         thickness = recHitTools_->getSiThickness(rh_detid);
2289       } else if (rh_detid.det() == DetId::HGCalHSc) {
2290         thickness = -1;
2291       } else {
2292         LogDebug("HGCalValidator") << "These are HGCal layer clusters, you shouldn't be here !!! " << layerid << "\n";
2293         continue;
2294       }
2295 
2296       //Count here only once the layer cluster and save the combination thick_layerid
2297       std::string curistr = std::to_string((int)thickness);
2298       std::string lay_string = std::to_string(layerid);
2299       while (lay_string.size() < 2)
2300         lay_string.insert(0, "0");
2301       curistr += "_" + lay_string;
2302       if (cluslay) {
2303         tnlcpl[layerid]++;
2304         istr = curistr;
2305         cluslay = false;
2306       }
2307 
2308       if ((thickness == 120.) && (recHitTools_->zside(rh_detid) > 0.)) {
2309         nthhits120p++;
2310       } else if ((thickness == 120.) && (recHitTools_->zside(rh_detid) < 0.)) {
2311         nthhits120m++;
2312       } else if ((thickness == 200.) && (recHitTools_->zside(rh_detid) > 0.)) {
2313         nthhits200p++;
2314       } else if ((thickness == 200.) && (recHitTools_->zside(rh_detid) < 0.)) {
2315         nthhits200m++;
2316       } else if ((thickness == 300.) && (recHitTools_->zside(rh_detid) > 0.)) {
2317         nthhits300p++;
2318       } else if ((thickness == 300.) && (recHitTools_->zside(rh_detid) < 0.)) {
2319         nthhits300m++;
2320       } else if ((thickness == -1) && (recHitTools_->zside(rh_detid) > 0.)) {
2321         nthhitsscintp++;
2322       } else if ((thickness == -1) && (recHitTools_->zside(rh_detid) < 0.)) {
2323         nthhitsscintm++;
2324       } else {  //assert(0);
2325         LogDebug("HGCalValidator")
2326             << " You are running a geometry that contains thicknesses different than the normal ones. "
2327             << "\n";
2328       }
2329 
2330       std::unordered_map<DetId, const HGCRecHit*>::const_iterator itcheck = hitMap.find(rh_detid);
2331       if (itcheck == hitMap.end()) {
2332         std::ostringstream st1;
2333         if ((rh_detid.det() == DetId::HGCalEE) || (rh_detid.det() == DetId::HGCalHSi)) {
2334           st1 << HGCSiliconDetId(rh_detid);
2335         } else if (rh_detid.det() == DetId::HGCalHSc) {
2336           st1 << HGCScintillatorDetId(rh_detid);
2337         } else {
2338           st1 << HFNoseDetId(rh_detid);
2339         }
2340         LogDebug("HGCalValidator") << " You shouldn't be here - Unable to find a hit " << rh_detid.rawId() << " "
2341                                    << rh_detid.det() << " " << st1.str() << "\n";
2342         continue;
2343       }
2344 
2345       const HGCRecHit* hit = itcheck->second;
2346 
2347       //Here for the per cell plots
2348       //----
2349       const double hit_x = recHitTools_->getPosition(rh_detid).x();
2350       const double hit_y = recHitTools_->getPosition(rh_detid).y();
2351       double distancetoseed = distance(seedx, seedy, hit_x, hit_y);
2352       double distancetomax = distance(maxx, maxy, hit_x, hit_y);
2353       if (distancetoseed != 0. && histograms.h_distancetoseedcell_perthickperlayer.count(curistr)) {
2354         histograms.h_distancetoseedcell_perthickperlayer.at(curistr)->Fill(distancetoseed);
2355       }
2356       //----
2357       if (distancetoseed != 0. && histograms.h_distancetoseedcell_perthickperlayer_eneweighted.count(curistr)) {
2358         histograms.h_distancetoseedcell_perthickperlayer_eneweighted.at(curistr)->Fill(distancetoseed, hit->energy());
2359       }
2360       //----
2361       if (distancetomax != 0. && histograms.h_distancetomaxcell_perthickperlayer.count(curistr)) {
2362         histograms.h_distancetomaxcell_perthickperlayer.at(curistr)->Fill(distancetomax);
2363       }
2364       //----
2365       if (distancetomax != 0. && histograms.h_distancetomaxcell_perthickperlayer_eneweighted.count(curistr)) {
2366         histograms.h_distancetomaxcell_perthickperlayer_eneweighted.at(curistr)->Fill(distancetomax, hit->energy());
2367       }
2368 
2369       //Let's check the density
2370       std::map<DetId, float>::const_iterator dit = densities.find(rh_detid);
2371       if (dit == densities.end()) {
2372         std::ostringstream st1;
2373         if ((rh_detid.det() == DetId::HGCalEE) || (rh_detid.det() == DetId::HGCalHSi)) {
2374           st1 << HGCSiliconDetId(rh_detid);
2375         } else if (rh_detid.det() == DetId::HGCalHSc) {
2376           st1 << HGCScintillatorDetId(rh_detid);
2377         } else {
2378           st1 << HFNoseDetId(rh_detid);
2379         }
2380         LogDebug("HGCalValidator") << " You shouldn't be here - Unable to find a density " << rh_detid.rawId() << " "
2381                                    << rh_detid.det() << " " << st1.str() << "\n";
2382         continue;
2383       }
2384 
2385       if (histograms.h_cellsenedens_perthick.count((int)thickness)) {
2386         histograms.h_cellsenedens_perthick.at((int)thickness)->Fill(dit->second);
2387       }
2388 
2389     }  // end of loop through hits and fractions
2390 
2391     //Check for simultaneously having hits of different kind. Checking at least two combinations is sufficient.
2392     if ((nthhits120p != 0 && nthhits200p != 0) || (nthhits120p != 0 && nthhits300p != 0) ||
2393         (nthhits120p != 0 && nthhitsscintp != 0) || (nthhits200p != 0 && nthhits300p != 0) ||
2394         (nthhits200p != 0 && nthhitsscintp != 0) || (nthhits300p != 0 && nthhitsscintp != 0)) {
2395       tnlcpthplus["mixed"]++;
2396     } else if ((nthhits120p != 0 || nthhits200p != 0 || nthhits300p != 0 || nthhitsscintp != 0)) {
2397       //This is a cluster with hits of one kind
2398       tnlcpthplus[std::to_string((int)thickness)]++;
2399     }
2400     if ((nthhits120m != 0 && nthhits200m != 0) || (nthhits120m != 0 && nthhits300m != 0) ||
2401         (nthhits120m != 0 && nthhitsscintm != 0) || (nthhits200m != 0 && nthhits300m != 0) ||
2402         (nthhits200m != 0 && nthhitsscintm != 0) || (nthhits300m != 0 && nthhitsscintm != 0)) {
2403       tnlcpthminus["mixed"]++;
2404     } else if ((nthhits120m != 0 || nthhits200m != 0 || nthhits300m != 0 || nthhitsscintm != 0)) {
2405       //This is a cluster with hits of one kind
2406       tnlcpthminus[std::to_string((int)thickness)]++;
2407     }
2408 
2409     //To find the thickness with the biggest amount of cells
2410     std::vector<int> bigamoth;
2411     bigamoth.clear();
2412     if (zside > 0) {
2413       bigamoth.push_back(nthhits120p);
2414       bigamoth.push_back(nthhits200p);
2415       bigamoth.push_back(nthhits300p);
2416       bigamoth.push_back(nthhitsscintp);
2417     } else if (zside < 0) {
2418       bigamoth.push_back(nthhits120m);
2419       bigamoth.push_back(nthhits200m);
2420       bigamoth.push_back(nthhits300m);
2421       bigamoth.push_back(nthhitsscintm);
2422     }
2423     auto bgth = std::max_element(bigamoth.begin(), bigamoth.end());
2424     istr = std::to_string(thicknesses[std::distance(bigamoth.begin(), bgth)]);
2425     std::string lay_string = std::to_string(layerid);
2426     while (lay_string.size() < 2)
2427       lay_string.insert(0, "0");
2428     istr += "_" + lay_string;
2429 
2430     //Here for the per cluster plots that need the thickness_layer info
2431     if (histograms.h_cellsnum_perthickperlayer.count(istr)) {
2432       histograms.h_cellsnum_perthickperlayer.at(istr)->Fill(hits_and_fractions.size());
2433     }
2434 
2435     //Now, with the distance between seed and max cell.
2436     double distancebetseedandmax = distance(seedx, seedy, maxx, maxy);
2437     //The thickness_layer combination in this case will use the thickness of the seed as a convention.
2438     std::string seedstr = std::to_string((int)recHitTools_->getSiThickness(seedid)) + "_" + std::to_string(layerid);
2439     seedstr += "_" + lay_string;
2440     if (histograms.h_distancebetseedandmaxcell_perthickperlayer.count(seedstr)) {
2441       histograms.h_distancebetseedandmaxcell_perthickperlayer.at(seedstr)->Fill(distancebetseedandmax);
2442     }
2443     const auto lc_en = lcId.energy();
2444     if (histograms.h_distancebetseedandmaxcellvsclusterenergy_perthickperlayer.count(seedstr)) {
2445       histograms.h_distancebetseedandmaxcellvsclusterenergy_perthickperlayer.at(seedstr)->Fill(distancebetseedandmax,
2446                                                                                                lc_en);
2447     }
2448 
2449     //Energy clustered per layer
2450     tecpl[layerid] = tecpl[layerid] + lc_en;
2451     ldbar[layerid] = ldbar[layerid] + lc_en * cummatbudg[(double)lay];
2452 
2453   }  //end of loop through clusters of the event
2454 
2455   // First a couple of variables to keep the sum of the energy of all clusters
2456   double sumeneallcluspl = 0.;
2457   double sumeneallclusmi = 0.;
2458   // and the longitudinal variable
2459   double sumldbarpl = 0.;
2460   double sumldbarmi = 0.;
2461   //Per layer : Loop 0->103
2462   for (unsigned ilayer = 0; ilayer < layers * 2; ++ilayer) {
2463     if (histograms.h_clusternum_perlayer.count(ilayer)) {
2464       histograms.h_clusternum_perlayer.at(ilayer)->Fill(tnlcpl[ilayer]);
2465     }
2466     // Two times one for plus and one for minus
2467     //First with the -z endcap
2468     if (ilayer < layers) {
2469       if (histograms.h_energyclustered_perlayer.count(ilayer)) {
2470         if (caloparteneminus != 0.) {
2471           histograms.h_energyclustered_perlayer.at(ilayer)->Fill(100. * tecpl[ilayer] / caloparteneminus);
2472         }
2473       }
2474       //Keep here the total energy for the event in -z
2475       sumeneallclusmi = sumeneallclusmi + tecpl[ilayer];
2476       //And for the longitudinal variable
2477       sumldbarmi = sumldbarmi + ldbar[ilayer];
2478     } else {  //Then for the +z
2479       if (histograms.h_energyclustered_perlayer.count(ilayer)) {
2480         if (caloparteneplus != 0.) {
2481           histograms.h_energyclustered_perlayer.at(ilayer)->Fill(100. * tecpl[ilayer] / caloparteneplus);
2482         }
2483       }
2484       //Keep here the total energy for the event in -z
2485       sumeneallcluspl = sumeneallcluspl + tecpl[ilayer];
2486       //And for the longitudinal variable
2487       sumldbarpl = sumldbarpl + ldbar[ilayer];
2488     }  //end of +z loop
2489 
2490   }  //end of loop over layers
2491 
2492   //Per thickness
2493   for (std::vector<int>::iterator it = thicknesses.begin(); it != thicknesses.end(); ++it) {
2494     if (histograms.h_clusternum_perthick.count(*it)) {
2495       histograms.h_clusternum_perthick.at(*it)->Fill(tnlcpthplus[std::to_string(*it)]);
2496       histograms.h_clusternum_perthick.at(*it)->Fill(tnlcpthminus[std::to_string(*it)]);
2497     }
2498   }
2499   //Mixed thickness clusters
2500   histograms.h_mixedhitscluster_zplus[count]->Fill(tnlcpthplus["mixed"]);
2501   histograms.h_mixedhitscluster_zminus[count]->Fill(tnlcpthminus["mixed"]);
2502 
2503   //Total energy clustered from all layer clusters (fraction)
2504   if (caloparteneplus != 0.) {
2505     histograms.h_energyclustered_zplus[count]->Fill(100. * sumeneallcluspl / caloparteneplus);
2506   }
2507   if (caloparteneminus != 0.) {
2508     histograms.h_energyclustered_zminus[count]->Fill(100. * sumeneallclusmi / caloparteneminus);
2509   }
2510 
2511   //For the longitudinal depth barycenter
2512   histograms.h_longdepthbarycentre_zplus[count]->Fill(sumldbarpl / sumeneallcluspl);
2513   histograms.h_longdepthbarycentre_zminus[count]->Fill(sumldbarmi / sumeneallclusmi);
2514 }
2515 
2516 void HGVHistoProducerAlgo::tracksters_to_SimTracksters(
2517     const Histograms& histograms,
2518     const int count,
2519     const ticl::TracksterCollection& tracksters,
2520     const reco::CaloClusterCollection& layerClusters,
2521     const ticl::TracksterCollection& simTSs,
2522     const validationType valType,
2523     const ticl::TracksterCollection& simTSs_fromCP,
2524     const std::map<unsigned int, std::vector<unsigned int>>& cpToSc_SimTrackstersMap,
2525     std::vector<SimCluster> const& sC,
2526     const edm::ProductID& cPHandle_id,
2527     std::vector<CaloParticle> const& cP,
2528     std::vector<size_t> const& cPIndices,
2529     std::vector<size_t> const& cPSelectedIndices,
2530     std::unordered_map<DetId, const HGCRecHit*> const& hitMap,
2531     unsigned int layers) const {
2532   const auto nTracksters = tracksters.size();
2533   const auto nSimTracksters = simTSs.size();
2534 
2535   std::unordered_map<DetId, std::vector<HGVHistoProducerAlgo::detIdInfoInCluster>> detIdSimTSId_Map;
2536   std::unordered_map<DetId, std::vector<HGVHistoProducerAlgo::detIdInfoInTrackster>> detIdToTracksterId_Map;
2537   std::vector<int> tracksters_FakeMerge(nTracksters, 0);
2538   std::vector<int> tracksters_PurityDuplicate(nTracksters, 0);
2539 
2540   // This vector contains the ids of the SimTracksters contributing with at least one hit to the Trackster and the reconstruction error
2541   //stsInTrackster[trackster][STSids]
2542   //Connects a Trackster with all related SimTracksters.
2543   std::vector<std::vector<std::pair<unsigned int, float>>> stsInTrackster;
2544   stsInTrackster.resize(nTracksters);
2545 
2546   // cPOnLayer[caloparticle]:
2547   //1. the sum of all rechits energy times fraction of the relevant simhit related to that calo particle.
2548   //2. the hits and fractions of that calo particle.
2549   //3. the layer clusters with matched rechit id.
2550   std::unordered_map<int, caloParticleOnLayer> cPOnLayer;
2551   std::unordered_map<int, std::vector<caloParticleOnLayer>> sCOnLayer;
2552   //Consider CaloParticles coming from the hard scatterer, excluding the PU contribution.
2553   for (const auto cpIndex : cPIndices) {
2554     cPOnLayer[cpIndex].caloParticleId = cpIndex;
2555     cPOnLayer[cpIndex].energy = 0.f;
2556     cPOnLayer[cpIndex].hits_and_fractions.clear();
2557     const auto nSC_inCP = sC.size();
2558     sCOnLayer[cpIndex].resize(nSC_inCP);
2559     for (unsigned int iSC = 0; iSC < nSC_inCP; iSC++) {
2560       sCOnLayer[cpIndex][iSC].caloParticleId = cpIndex;
2561       sCOnLayer[cpIndex][iSC].energy = 0.f;
2562       sCOnLayer[cpIndex][iSC].hits_and_fractions.clear();
2563     }
2564   }
2565 
2566   auto getCPId = [](const ticl::Trackster& simTS,
2567                     const unsigned int iSTS,
2568                     const edm::ProductID& cPHandle_id,
2569                     const std::map<unsigned int, std::vector<unsigned int>>& cpToSc_SimTrackstersMap,
2570                     const ticl::TracksterCollection& simTSs_fromCP) {
2571     unsigned int cpId = -1;
2572 
2573     const auto productID = simTS.seedID();
2574     if (productID == cPHandle_id) {
2575       cpId = simTS.seedIndex();
2576     } else {  // SimTrackster from SimCluster
2577       const auto findSimTSFromCP = std::find_if(
2578           std::begin(cpToSc_SimTrackstersMap),
2579           std::end(cpToSc_SimTrackstersMap),
2580           [&](const std::pair<unsigned int, std::vector<unsigned int>>& cpToScs) {
2581             return std::find(std::begin(cpToScs.second), std::end(cpToScs.second), iSTS) != std::end(cpToScs.second);
2582           });
2583       if (findSimTSFromCP != std::end(cpToSc_SimTrackstersMap)) {
2584         cpId = simTSs_fromCP[findSimTSFromCP->first].seedIndex();
2585       }
2586     }
2587 
2588     return cpId;
2589   };
2590 
2591   auto getLCId = [](const std::vector<unsigned int>& tst_vertices,
2592                     const reco::CaloClusterCollection& layerClusters,
2593                     const DetId& hitid) {
2594     unsigned int lcId = -1;
2595     std::for_each(std::begin(tst_vertices), std::end(tst_vertices), [&](unsigned int idx) {
2596       const auto& lc_haf = layerClusters[idx].hitsAndFractions();
2597       const auto& hitFound = std::find_if(std::begin(lc_haf),
2598                                           std::end(lc_haf),
2599                                           [&hitid](const std::pair<DetId, float>& v) { return v.first == hitid; });
2600       if (hitFound != lc_haf.end())  // not all hits may be clusterized
2601         lcId = idx;
2602     });
2603     return lcId;
2604   };
2605 
2606   for (unsigned int iSTS = 0; iSTS < nSimTracksters; ++iSTS) {
2607     const auto cpId = getCPId(simTSs[iSTS], iSTS, cPHandle_id, cpToSc_SimTrackstersMap, simTSs_fromCP);
2608     if (std::find(cPIndices.begin(), cPIndices.end(), cpId) == cPIndices.end())
2609       continue;
2610 
2611     // Loop through SimClusters
2612     for (const auto& simCluster : cP[cpId].simClusters()) {
2613       auto iSim = simTSs[iSTS].seedIndex();
2614       if (simTSs[iSTS].seedID() != cPHandle_id) {  // SimTrackster from SimCluster
2615         if (iSim != (&(*simCluster) - &(sC[0])))
2616           continue;
2617       } else
2618         iSim = 0;
2619 
2620       for (const auto& it_haf : simCluster->hits_and_fractions()) {
2621         const auto hitid = (it_haf.first);
2622         const auto lcId = getLCId(simTSs[iSTS].vertices(), layerClusters, hitid);
2623         //V9:maps the layers in -z: 0->51 and in +z: 52->103
2624         //V10:maps the layers in -z: 0->49 and in +z: 50->99
2625         const auto itcheck = hitMap.find(hitid);
2626         //Check whether the current hit belonging to sim cluster has a reconstructed hit
2627         if ((valType == 0 && itcheck != hitMap.end()) || (valType > 0 && int(lcId) >= 0)) {
2628           float lcFraction = 0;
2629           if (valType > 0) {
2630             const auto iLC = std::find(simTSs[iSTS].vertices().begin(), simTSs[iSTS].vertices().end(), lcId);
2631             lcFraction =
2632                 1.f / simTSs[iSTS].vertex_multiplicity(std::distance(std::begin(simTSs[iSTS].vertices()), iLC));
2633           }
2634           const auto elemId = (valType == 0) ? hitid : lcId;
2635           const auto elemFr = (valType == 0) ? it_haf.second : lcFraction;
2636           //Since the current hit from SimCluster has a reconstructed hit {with the same detid, belonging to the corresponding SimTrackster},
2637           //make a map that will connect a {detid,lcID} with:
2638           //1. the SimTracksters that have a SimCluster with {sim hits in, LCs containing} that cell via SimTrackster id.
2639           //2. the sum of all {SimHits, LCs} fractions that contributes to that detid.
2640           //So, keep in mind that in case of multiple CaloParticles contributing in the same cell
2641           //the fraction is the sum over all calo particles. So, something like:
2642           //detid: (caloparticle 1, sum of hits fractions in that detid over all cp) , (caloparticle 2, sum of hits fractions in that detid over all cp), (caloparticle 3, sum of hits fractions in that detid over all cp) ...
2643           if (detIdSimTSId_Map.find(elemId) == detIdSimTSId_Map.end()) {
2644             detIdSimTSId_Map[elemId] = std::vector<HGVHistoProducerAlgo::detIdInfoInCluster>();
2645             detIdSimTSId_Map[elemId].emplace_back(HGVHistoProducerAlgo::detIdInfoInCluster{iSTS, elemFr});
2646           } else {
2647             auto findSTSIt =
2648                 std::find(detIdSimTSId_Map[elemId].begin(),
2649                           detIdSimTSId_Map[elemId].end(),
2650                           HGVHistoProducerAlgo::detIdInfoInCluster{
2651                               iSTS, 0});  // only the first element is used for the matching (overloaded operator==)
2652             if (findSTSIt != detIdSimTSId_Map[elemId].end()) {
2653               if (valType == 0)
2654                 findSTSIt->fraction += elemFr;
2655             } else {
2656               detIdSimTSId_Map[elemId].emplace_back(HGVHistoProducerAlgo::detIdInfoInCluster{iSTS, elemFr});
2657             }
2658           }
2659           const auto hitEn = itcheck->second->energy();
2660           //Since the current hit from sim cluster has a reconstructed hit with the same detid,
2661           //fill the cPOnLayer[caloparticle][layer] object with energy (sum of all rechits energy times fraction
2662           //of the relevant simhit) and keep the hit (detid and fraction) that contributed.
2663           cPOnLayer[cpId].energy += it_haf.second * hitEn;
2664           sCOnLayer[cpId][iSim].energy += elemFr * hitEn;
2665           // Need to compress the hits and fractions in order to have a
2666           // reasonable score between CP and LC. Imagine, for example, that a
2667           // CP has detID X used by 2 SimClusters with different fractions. If
2668           // a single LC uses X with fraction 1 and is compared to the 2
2669           // contributions separately, it will be assigned a score != 0, which
2670           // is wrong.
2671           auto& haf = cPOnLayer[cpId].hits_and_fractions;
2672           auto found = std::find_if(
2673               std::begin(haf), std::end(haf), [&hitid](const std::pair<DetId, float>& v) { return v.first == hitid; });
2674           if (found != haf.end())
2675             found->second += it_haf.second;
2676           else
2677             haf.emplace_back(hitid, it_haf.second);
2678           // Same for sCOnLayer
2679           auto& haf_sc = sCOnLayer[cpId][iSim].hits_and_fractions;
2680           auto found_sc = std::find_if(std::begin(haf_sc),
2681                                        std::end(haf_sc),
2682                                        [&hitid](const std::pair<DetId, float>& v) { return v.first == hitid; });
2683           if (found_sc != haf_sc.end())
2684             found_sc->second += it_haf.second;
2685           else
2686             haf_sc.emplace_back(hitid, it_haf.second);
2687         }
2688       }  // end of loop through SimHits
2689     }    // end of loop through SimClusters
2690   }      // end of loop through SimTracksters
2691 
2692   auto apply_LCMultiplicity = [](const ticl::Trackster& trackster, const reco::CaloClusterCollection& layerClusters) {
2693     std::vector<std::pair<DetId, float>> hits_and_fractions_norm;
2694     int lcInTst = 0;
2695     std::for_each(std::begin(trackster.vertices()), std::end(trackster.vertices()), [&](unsigned int idx) {
2696       const auto fraction = 1.f / trackster.vertex_multiplicity(lcInTst++);
2697       for (const auto& cell : layerClusters[idx].hitsAndFractions()) {
2698         hits_and_fractions_norm.emplace_back(
2699             cell.first, cell.second * fraction);  // cell.second is the hit fraction in the layerCluster
2700       }
2701     });
2702     return hits_and_fractions_norm;
2703   };
2704 
2705   auto ScoreCutSTStoTSPurDup = ScoreCutSTStoTSPurDup_[0];
2706   auto ScoreCutTStoSTSFakeMerge = ScoreCutTStoSTSFakeMerge_[0];
2707   // Loop through Tracksters
2708   for (unsigned int tstId = 0; tstId < nTracksters; ++tstId) {
2709     const auto& tst = tracksters[tstId];
2710     if (tstId == 0)
2711       if ((valType > 0) && (tst.ticlIteration() == ticl::Trackster::SIM)) {
2712         ScoreCutSTStoTSPurDup = ScoreCutSTStoTSPurDup_[valType];
2713         ScoreCutTStoSTSFakeMerge = ScoreCutTStoSTSFakeMerge_[valType];
2714       }
2715 
2716     if (tst.vertices().empty())
2717       continue;
2718 
2719     std::unordered_map<unsigned, float> CPEnergyInTS;
2720     int maxCPId_byNumberOfHits = -1;
2721     unsigned int maxCPNumberOfHitsInTS = 0;
2722     int maxCPId_byEnergy = -1;
2723     float maxEnergySharedTSandCP = 0.f;
2724     float energyFractionOfTSinCP = 0.f;
2725     float energyFractionOfCPinTS = 0.f;
2726 
2727     //In case of matched rechit-simhit, so matched
2728     //CaloParticle-LayerCluster-Trackster, he counts and saves the number of
2729     //rechits related to the maximum energy CaloParticle out of all
2730     //CaloParticles related to that layer cluster and Trackster.
2731 
2732     std::unordered_map<unsigned, unsigned> occurrencesCPinTS;
2733     unsigned int numberOfNoiseHitsInTS = 0;
2734     unsigned int numberOfHaloHitsInTS = 0;
2735 
2736     const auto tst_hitsAndFractions = apply_LCMultiplicity(tst, layerClusters);
2737     const auto numberOfHitsInTS = tst_hitsAndFractions.size();
2738 
2739     //hitsToCaloParticleId is a vector of ints, one for each rechit of the
2740     //layer cluster under study. If negative, there is no simhit from any CaloParticle related.
2741     //If positive, at least one CaloParticle has been found with matched simhit.
2742     //In more detail:
2743     // 1. hitsToCaloParticleId[iHit] = -3
2744     //    TN:  These represent Halo Cells(N) that have not been
2745     //    assigned to any CaloParticle (hence the T).
2746     // 2. hitsToCaloParticleId[iHit] = -2
2747     //    FN: There represent Halo Cells(N) that have been assigned
2748     //    to a CaloParticle (hence the F, since those should have not been marked as halo)
2749     // 3. hitsToCaloParticleId[iHit] = -1
2750     //    FP: These represent Real Cells(P) that have not been
2751     //    assigned to any CaloParticle (hence the F, since these are fakes)
2752     // 4. hitsToCaloParticleId[iHit] >= 0
2753     //    TP There represent Real Cells(P) that have been assigned
2754     //    to a CaloParticle (hence the T)
2755     std::vector<int> hitsToCaloParticleId(numberOfHitsInTS);
2756 
2757     //Loop through the hits of the trackster under study
2758     for (unsigned int iHit = 0; iHit < numberOfHitsInTS; iHit++) {
2759       const auto rh_detid = tst_hitsAndFractions[iHit].first;
2760       const auto rhFraction = tst_hitsAndFractions[iHit].second;
2761 
2762       const auto lcId_r = getLCId(tst.vertices(), layerClusters, rh_detid);
2763       const auto iLC_r = std::find(tst.vertices().begin(), tst.vertices().end(), lcId_r);
2764       const auto lcFraction_r = 1.f / tst.vertex_multiplicity(std::distance(std::begin(tst.vertices()), iLC_r));
2765 
2766       //Make a map that will connect a detid (that belongs to a rechit of the layer cluster under study,
2767       //no need to save others) with:
2768       //1. the layer clusters that have rechits in that detid
2769       //2. the fraction of the rechit of each layer cluster that contributes to that detid.
2770       //So, something like:
2771       //detid: (layer cluster 1, hit fraction) , (layer cluster 2, hit fraction), (layer cluster 3, hit fraction) ...
2772       //here comparing with the calo particle map above the
2773       if (detIdToTracksterId_Map.find(rh_detid) == detIdToTracksterId_Map.end()) {
2774         detIdToTracksterId_Map[rh_detid] = std::vector<HGVHistoProducerAlgo::detIdInfoInTrackster>();
2775         detIdToTracksterId_Map[rh_detid].emplace_back(
2776             HGVHistoProducerAlgo::detIdInfoInTrackster{tstId, lcId_r, rhFraction});
2777       } else {
2778         auto findTSIt =
2779             std::find(detIdToTracksterId_Map[rh_detid].begin(),
2780                       detIdToTracksterId_Map[rh_detid].end(),
2781                       HGVHistoProducerAlgo::detIdInfoInTrackster{
2782                           tstId, 0, 0});  // only the first element is used for the matching (overloaded operator==)
2783         if (findTSIt != detIdToTracksterId_Map[rh_detid].end()) {
2784           if (valType == 0)
2785             findTSIt->fraction += rhFraction;
2786         } else {
2787           detIdToTracksterId_Map[rh_detid].emplace_back(
2788               HGVHistoProducerAlgo::detIdInfoInTrackster{tstId, lcId_r, rhFraction});
2789         }
2790       }
2791 
2792       // if the fraction is zero or the hit does not belong to any calo
2793       // particle, set the caloparticleId for the hit to -1 this will
2794       // contribute to the number of noise hits
2795       // MR Remove the case in which the fraction is 0, since this could be a
2796       // real hit that has been marked as halo.
2797       if (rhFraction == 0.) {
2798         hitsToCaloParticleId[iHit] = -2;
2799         numberOfHaloHitsInTS++;
2800       }
2801 
2802       // Check whether the RecHit of the trackster under study has a SimHit in the same cell
2803       const auto elemId = (valType == 0) ? rh_detid.rawId() : lcId_r;
2804       const auto recoFr = (valType == 0) ? rhFraction : lcFraction_r;
2805       const auto& hit_find_in_STS = detIdSimTSId_Map.find(elemId);
2806       if (hit_find_in_STS == detIdSimTSId_Map.end()) {
2807         hitsToCaloParticleId[iHit] -= 1;
2808       } else {
2809         // Since the hit is belonging to the layer cluster, it must be also in the rechits map
2810         const auto hitEn = hitMap.find(rh_detid)->second->energy();
2811         //const auto layerId =
2812         //recHitTools_->getLayerWithOffset(rh_detid) + layers * ((recHitTools_->zside(rh_detid) + 1) >> 1) - 1;
2813         //0;
2814 
2815         auto maxCPEnergyInTS = 0.f;
2816         auto maxCPId = -1;
2817         for (const auto& h : hit_find_in_STS->second) {
2818           const auto shared_fraction = std::min(recoFr, h.fraction);
2819           const auto iSTS = h.clusterId;
2820           const auto& simTS = simTSs[iSTS];
2821           auto iSim = simTS.seedIndex();
2822           if (simTSs[iSTS].seedID() == cPHandle_id)  // SimTrackster from CaloParticle
2823             iSim = 0;
2824 
2825           // SimTrackster with simHits connected via detid with the rechit under study
2826           //So, from all layers clusters, find the rechits that are connected with a calo particle and save/calculate the
2827           //energy of that calo particle as the sum over all rechits of the rechits energy weighted
2828           //by the caloparticle's fraction related to that rechit.
2829           const auto cpId = getCPId(simTS, iSTS, cPHandle_id, cpToSc_SimTrackstersMap, simTSs_fromCP);
2830           if (std::find(cPIndices.begin(), cPIndices.end(), cpId) == cPIndices.end())
2831             continue;
2832 
2833           CPEnergyInTS[cpId] += shared_fraction * hitEn;
2834           //Here cPOnLayer[caloparticle][layer] describe above is set.
2835           //Here for Tracksters with matched rechit the CP fraction times hit energy is added and saved .
2836           cPOnLayer[cpId].layerClusterIdToEnergyAndScore[tstId].first += shared_fraction * hitEn;
2837           sCOnLayer[cpId][iSim].layerClusterIdToEnergyAndScore[tstId].first += shared_fraction * hitEn;
2838           cPOnLayer[cpId].layerClusterIdToEnergyAndScore[tstId].second = FLT_MAX;
2839           sCOnLayer[cpId][iSim].layerClusterIdToEnergyAndScore[tstId].second = FLT_MAX;
2840           //stsInTrackster[trackster][STSids]
2841           //Connects a Trackster with all related SimTracksters.
2842           stsInTrackster[tstId].emplace_back(iSTS, FLT_MAX);
2843           //From all CaloParticles related to a layer cluster, it saves id and energy of the calo particle
2844           //that after simhit-rechit matching in layer has the maximum energy.
2845           if (shared_fraction > maxCPEnergyInTS) {
2846             //energy is used only here. cpid is saved for Tracksters
2847             maxCPEnergyInTS = CPEnergyInTS[cpId];
2848             maxCPId = cpId;
2849           }
2850         }
2851         //Keep in mind here maxCPId could be zero. So, below ask for negative not including zero to count noise.
2852         hitsToCaloParticleId[iHit] = maxCPId;
2853       }
2854 
2855     }  //end of loop through rechits of the layer cluster.
2856 
2857     //Loop through all rechits to count how many of them are noise and how many are matched.
2858     //In case of matched rechit-simhit, he counts and saves the number of rechits related to the maximum energy CaloParticle.
2859     for (auto c : hitsToCaloParticleId) {
2860       if (c < 0)
2861         numberOfNoiseHitsInTS++;
2862       else
2863         occurrencesCPinTS[c]++;
2864     }
2865 
2866     //Below from all maximum energy CaloParticles, he saves the one with the largest amount
2867     //of related rechits.
2868     for (auto& c : occurrencesCPinTS) {
2869       if (c.second > maxCPNumberOfHitsInTS) {
2870         maxCPId_byNumberOfHits = c.first;
2871         maxCPNumberOfHitsInTS = c.second;
2872       }
2873     }
2874 
2875     //Find the CaloParticle that has the maximum energy shared with the Trackster under study.
2876     for (auto& c : CPEnergyInTS) {
2877       if (c.second > maxEnergySharedTSandCP) {
2878         maxCPId_byEnergy = c.first;
2879         maxEnergySharedTSandCP = c.second;
2880       }
2881     }
2882     //The energy of the CaloParticle that found to have the maximum energy shared with the Trackster under study.
2883     float totalCPEnergyFromLayerCP = 0.f;
2884     if (maxCPId_byEnergy >= 0) {
2885       totalCPEnergyFromLayerCP += cPOnLayer[maxCPId_byEnergy].energy;
2886       energyFractionOfCPinTS = maxEnergySharedTSandCP / totalCPEnergyFromLayerCP;
2887       if (tst.raw_energy() > 0.f) {
2888         energyFractionOfTSinCP = maxEnergySharedTSandCP / tst.raw_energy();
2889       }
2890     }
2891 
2892     LogDebug("HGCalValidator") << std::setw(12) << "Trackster\t" << std::setw(10) << "energy\t" << std::setw(5)
2893                                << "nhits\t" << std::setw(12) << "noise hits\t" << std::setw(22)
2894                                << "maxCPId_byNumberOfHits\t" << std::setw(8) << "nhitsCP\t" << std::setw(16)
2895                                << "maxCPId_byEnergy\t" << std::setw(23) << "maxEnergySharedTSandCP\t" << std::setw(22)
2896                                << "totalCPEnergyFromAllLayerCP\t" << std::setw(22) << "energyFractionOfTSinCP\t"
2897                                << std::setw(25) << "energyFractionOfCPinTS\t" << std::endl;
2898     LogDebug("HGCalValidator") << std::setw(12) << tstId << "\t"  //LogDebug("HGCalValidator")
2899                                << std::setw(10) << tst.raw_energy() << "\t" << std::setw(5) << numberOfHitsInTS << "\t"
2900                                << std::setw(12) << numberOfNoiseHitsInTS << "\t" << std::setw(22)
2901                                << maxCPId_byNumberOfHits << "\t" << std::setw(8) << maxCPNumberOfHitsInTS << "\t"
2902                                << std::setw(16) << maxCPId_byEnergy << "\t" << std::setw(23) << maxEnergySharedTSandCP
2903                                << "\t" << std::setw(22) << totalCPEnergyFromLayerCP << "\t" << std::setw(22)
2904                                << energyFractionOfTSinCP << "\t" << std::setw(25) << energyFractionOfCPinTS
2905                                << std::endl;
2906 
2907   }  // end of loop through Tracksters
2908 
2909   // Loop through Tracksters
2910   for (unsigned int tstId = 0; tstId < nTracksters; ++tstId) {
2911     const auto& tst = tracksters[tstId];
2912     if (tst.vertices().empty())
2913       continue;
2914 
2915     // find the unique SimTrackster ids contributing to the Trackster
2916     //stsInTrackster[trackster][STSids]
2917     std::sort(stsInTrackster[tstId].begin(), stsInTrackster[tstId].end());
2918     const auto last = std::unique(stsInTrackster[tstId].begin(), stsInTrackster[tstId].end());
2919     stsInTrackster[tstId].erase(last, stsInTrackster[tstId].end());
2920 
2921     if (tst.raw_energy() == 0. && !stsInTrackster[tstId].empty()) {
2922       //Loop through all SimTracksters contributing to Trackster tstId
2923       for (auto& stsPair : stsInTrackster[tstId]) {
2924         // In case of a Trackster with zero energy but related SimTracksters the score is set to 1
2925         stsPair.second = 1.;
2926         LogDebug("HGCalValidator") << "Trackster Id:\t" << tstId << "\tSimTrackster id:\t" << stsPair.first
2927                                    << "\tscore\t" << stsPair.second << std::endl;
2928         histograms.h_score_trackster2caloparticle[valType][count]->Fill(stsPair.second);
2929       }
2930       continue;
2931     }
2932 
2933     const auto tst_hitsAndFractions = apply_LCMultiplicity(tst, layerClusters);
2934 
2935     // Compute the correct normalization
2936     float tracksterEnergy = 0.f, invTracksterEnergyWeight = 0.f;
2937     for (const auto& haf : tst_hitsAndFractions) {
2938       float hitFr = 0.f;
2939       if (valType == 0) {
2940         hitFr = haf.second;
2941       } else {
2942         const auto lcId = getLCId(tst.vertices(), layerClusters, haf.first);
2943         const auto iLC = std::find(tst.vertices().begin(), tst.vertices().end(), lcId);
2944         hitFr = 1.f / tst.vertex_multiplicity(std::distance(std::begin(tst.vertices()), iLC));
2945       }
2946       tracksterEnergy += hitFr * hitMap.at(haf.first)->energy();
2947       invTracksterEnergyWeight += pow(hitFr * hitMap.at(haf.first)->energy(), 2);
2948     }
2949     if (invTracksterEnergyWeight)
2950       invTracksterEnergyWeight = 1.f / invTracksterEnergyWeight;
2951 
2952     for (const auto& haf : tst_hitsAndFractions) {
2953       const auto rh_detid = haf.first;
2954       unsigned int elemId = 0;
2955       float rhFraction = 0.f;
2956       if (valType == 0) {
2957         elemId = rh_detid.rawId();
2958         rhFraction = haf.second;
2959       } else {
2960         const auto lcId = getLCId(tst.vertices(), layerClusters, rh_detid);
2961         elemId = lcId;
2962         const auto iLC = std::find(tst.vertices().begin(), tst.vertices().end(), lcId);
2963         rhFraction = 1.f / tst.vertex_multiplicity(std::distance(std::begin(tst.vertices()), iLC));
2964       }
2965 
2966       bool hitWithNoSTS = false;
2967       if (detIdSimTSId_Map.find(elemId) == detIdSimTSId_Map.end())
2968         hitWithNoSTS = true;
2969       const HGCRecHit* hit = hitMap.find(rh_detid)->second;
2970       const auto hitEnergyWeight = pow(hit->energy(), 2);
2971 
2972       for (auto& stsPair : stsInTrackster[tstId]) {
2973         float cpFraction = 0.f;
2974         if (!hitWithNoSTS) {
2975           const auto& findSTSIt = std::find(
2976               detIdSimTSId_Map[elemId].begin(),
2977               detIdSimTSId_Map[elemId].end(),
2978               HGVHistoProducerAlgo::detIdInfoInCluster{
2979                   stsPair.first, 0.f});  // only the first element is used for the matching (overloaded operator==)
2980           if (findSTSIt != detIdSimTSId_Map[elemId].end())
2981             cpFraction = findSTSIt->fraction;
2982         }
2983         if (stsPair.second == FLT_MAX) {
2984           stsPair.second = 0.f;
2985         }
2986         stsPair.second +=
2987             min(pow(rhFraction - cpFraction, 2), pow(rhFraction, 2)) * hitEnergyWeight * invTracksterEnergyWeight;
2988       }
2989     }  // end of loop through trackster rechits
2990 
2991     //In case of a Trackster with some energy but none related CaloParticles print some info.
2992     if (stsInTrackster[tstId].empty())
2993       LogDebug("HGCalValidator") << "Trackster Id: " << tstId << "\tSimTrackster id: -1"
2994                                  << "\tscore: -1\n";
2995 
2996     tracksters_FakeMerge[tstId] =
2997         std::count_if(std::begin(stsInTrackster[tstId]),
2998                       std::end(stsInTrackster[tstId]),
2999                       [ScoreCutTStoSTSFakeMerge](const auto& obj) { return obj.second < ScoreCutTStoSTSFakeMerge; });
3000 
3001     const auto score = std::min_element(std::begin(stsInTrackster[tstId]),
3002                                         std::end(stsInTrackster[tstId]),
3003                                         [](const auto& obj1, const auto& obj2) { return obj1.second < obj2.second; });
3004     float score2 = -1;
3005     float sharedEneFrac2 = 0;
3006     for (const auto& stsPair : stsInTrackster[tstId]) {
3007       const auto iSTS = stsPair.first;
3008       const auto iScore = stsPair.second;
3009       const auto cpId = getCPId(simTSs[iSTS], iSTS, cPHandle_id, cpToSc_SimTrackstersMap, simTSs_fromCP);
3010       auto iSim = simTSs[iSTS].seedIndex();
3011       if (simTSs[iSTS].seedID() == cPHandle_id)  // SimTrackster from CaloParticle
3012         iSim = 0;
3013       const auto& simOnLayer = (valType == 0) ? cPOnLayer[cpId] : sCOnLayer[cpId][iSim];
3014 
3015       float sharedeneCPallLayers = 0.;
3016       sharedeneCPallLayers += simOnLayer.layerClusterIdToEnergyAndScore.count(tstId)
3017                                   ? simOnLayer.layerClusterIdToEnergyAndScore.at(tstId).first
3018                                   : 0;
3019       if (tracksterEnergy == 0)
3020         continue;
3021       const auto sharedEneFrac = sharedeneCPallLayers / tracksterEnergy;
3022       LogDebug("HGCalValidator") << "\nTrackster id: " << tstId << " (" << tst.vertices().size() << " vertices)"
3023                                  << "\tSimTrackster Id: " << iSTS << " (" << simTSs[iSTS].vertices().size()
3024                                  << " vertices)"
3025                                  << " (CP id: " << cpId << ")\tscore: " << iScore
3026                                  << "\tsharedeneCPallLayers: " << sharedeneCPallLayers << std::endl;
3027 
3028       histograms.h_score_trackster2caloparticle[valType][count]->Fill(iScore);
3029       histograms.h_sharedenergy_trackster2caloparticle[valType][count]->Fill(sharedEneFrac);
3030       histograms.h_energy_vs_score_trackster2caloparticle[valType][count]->Fill(iScore, sharedEneFrac);
3031       if (iSTS == score->first) {
3032         histograms.h_score_trackster2bestCaloparticle[valType][count]->Fill(iScore);
3033         histograms.h_sharedenergy_trackster2bestCaloparticle[valType][count]->Fill(sharedEneFrac);
3034         histograms.h_sharedenergy_trackster2bestCaloparticle_vs_eta[valType][count]->Fill(tst.barycenter().eta(),
3035                                                                                           sharedEneFrac);
3036         histograms.h_sharedenergy_trackster2bestCaloparticle_vs_phi[valType][count]->Fill(tst.barycenter().phi(),
3037                                                                                           sharedEneFrac);
3038         histograms.h_energy_vs_score_trackster2bestCaloparticle[valType][count]->Fill(iScore, sharedEneFrac);
3039       } else if (score2 < 0 || iScore < score2) {
3040         score2 = iScore;
3041         sharedEneFrac2 = sharedEneFrac;
3042       }
3043     }  // end of loop through SimTracksters associated to Trackster
3044     if (score2 > -1) {
3045       histograms.h_score_trackster2bestCaloparticle2[valType][count]->Fill(score2);
3046       histograms.h_sharedenergy_trackster2bestCaloparticle2[valType][count]->Fill(sharedEneFrac2);
3047       histograms.h_energy_vs_score_trackster2bestCaloparticle2[valType][count]->Fill(score2, sharedEneFrac2);
3048     }
3049   }  // end of loop through Tracksters
3050 
3051   std::unordered_map<unsigned int, std::vector<float>> score3d;
3052   std::unordered_map<unsigned int, std::vector<float>> tstSharedEnergy;
3053 
3054   for (unsigned int iSTS = 0; iSTS < nSimTracksters; ++iSTS) {
3055     score3d[iSTS].resize(nTracksters);
3056     tstSharedEnergy[iSTS].resize(nTracksters);
3057     for (unsigned int j = 0; j < nTracksters; ++j) {
3058       score3d[iSTS][j] = FLT_MAX;
3059       tstSharedEnergy[iSTS][j] = 0.f;
3060     }
3061   }
3062 
3063   // Fill the plots to compute the different metrics linked to
3064   // gen-level, namely efficiency, purity and duplicate. In this loop should restrict
3065   // only to the selected caloParaticles.
3066   for (unsigned int iSTS = 0; iSTS < nSimTracksters; ++iSTS) {
3067     const auto& sts = simTSs[iSTS];
3068     const auto& cpId = getCPId(sts, iSTS, cPHandle_id, cpToSc_SimTrackstersMap, simTSs_fromCP);
3069     if (valType == 0 && std::find(cPSelectedIndices.begin(), cPSelectedIndices.end(), cpId) == cPSelectedIndices.end())
3070       continue;
3071 
3072     const auto& hafLC = apply_LCMultiplicity(sts, layerClusters);
3073     float SimEnergy_LC = 0.f;
3074     for (const auto& haf : hafLC) {
3075       const auto lcId = getLCId(sts.vertices(), layerClusters, haf.first);
3076       const auto iLC = std::find(sts.vertices().begin(), sts.vertices().end(), lcId);
3077       SimEnergy_LC +=
3078           hitMap.at(haf.first)->energy() / sts.vertex_multiplicity(std::distance(std::begin(sts.vertices()), iLC));
3079     }
3080 
3081     auto iSim = sts.seedIndex();
3082     if (sts.seedID() == cPHandle_id)  // SimTrackster from CaloParticle
3083       iSim = 0;
3084     auto& simOnLayer = (valType == 0) ? cPOnLayer[cpId] : sCOnLayer[cpId][iSim];
3085 
3086     // Keep the Trackster ids that are related to
3087     // SimTrackster under study for the final filling of the score
3088     std::set<unsigned int> stsId_tstId_related;
3089     auto& score3d_iSTS = score3d[iSTS];
3090 
3091     float SimEnergy = 0.f;
3092     float SimEnergyWeight = 0.f, hitsEnergyWeight = 0.f;
3093     //for (unsigned int layerId = 0; layerId < 1/*layers * 2*/; ++layerId) {
3094     const auto SimNumberOfHits = simOnLayer.hits_and_fractions.size();
3095     if (SimNumberOfHits == 0)
3096       continue;
3097     SimEnergy += simOnLayer.energy;
3098     int tstWithMaxEnergyInCP = -1;
3099     //This is the maximum energy related to Trackster per layer.
3100     float maxEnergyTSperlayerinSim = 0.f;
3101     float SimEnergyFractionInTSperlayer = 0.f;
3102     //Remember and not confused by name. layerClusterIdToEnergyAndScore contains the Trackster id.
3103     for (const auto& tst : simOnLayer.layerClusterIdToEnergyAndScore) {
3104       if (tst.second.first > maxEnergyTSperlayerinSim) {
3105         maxEnergyTSperlayerinSim = tst.second.first;
3106         tstWithMaxEnergyInCP = tst.first;
3107       }
3108     }
3109     if (SimEnergy > 0.f)
3110       SimEnergyFractionInTSperlayer = maxEnergyTSperlayerinSim / SimEnergy;
3111 
3112     LogDebug("HGCalValidator") << std::setw(12) << "caloparticle\t" << std::setw(15) << "cp total energy\t"
3113                                << std::setw(15) << "cpEnergyOnLayer\t" << std::setw(14) << "CPNhitsOnLayer\t"
3114                                << std::setw(18) << "tstWithMaxEnergyInCP\t" << std::setw(15) << "maxEnergyTSinCP\t"
3115                                << std::setw(20) << "CPEnergyFractionInTS"
3116                                << "\n";
3117     LogDebug("HGCalValidator") << std::setw(12) << cpId << "\t" << std::setw(15) << sts.raw_energy() << "\t"
3118                                << std::setw(15) << SimEnergy << "\t" << std::setw(14) << SimNumberOfHits << "\t"
3119                                << std::setw(18) << tstWithMaxEnergyInCP << "\t" << std::setw(15)
3120                                << maxEnergyTSperlayerinSim << "\t" << std::setw(20) << SimEnergyFractionInTSperlayer
3121                                << "\n";
3122 
3123     for (const auto& haf : ((valType == 0) ? simOnLayer.hits_and_fractions : hafLC)) {
3124       const auto& hitDetId = haf.first;
3125       // Compute the correct normalization
3126       // Need to loop on the simOnLayer data structure since this is the
3127       // only one that has the compressed information for multiple usage
3128       // of the same DetId by different SimClusters by a single CaloParticle.
3129       SimEnergyWeight += pow(haf.second * hitMap.at(hitDetId)->energy(), 2);
3130 
3131       const auto lcId = getLCId(sts.vertices(), layerClusters, hitDetId);
3132       float cpFraction = 0.f;
3133       if (valType == 0) {
3134         cpFraction = haf.second;
3135       } else {
3136         const auto iLC = std::find(sts.vertices().begin(), sts.vertices().end(), lcId);
3137         cpFraction = 1.f / sts.vertex_multiplicity(std::distance(std::begin(sts.vertices()), iLC));
3138       }
3139       if (cpFraction == 0.f)
3140         continue;  // hopefully this should never happen
3141 
3142       bool hitWithNoTS = false;
3143       if (detIdToTracksterId_Map.find(hitDetId) == detIdToTracksterId_Map.end())
3144         hitWithNoTS = true;
3145       const HGCRecHit* hit = hitMap.find(hitDetId)->second;
3146       const auto hitEnergyWeight = pow(hit->energy(), 2);
3147       hitsEnergyWeight += pow(cpFraction, 2) * hitEnergyWeight;
3148 
3149       for (auto& tsPair : simOnLayer.layerClusterIdToEnergyAndScore) {
3150         const auto tstId = tsPair.first;
3151         stsId_tstId_related.insert(tstId);
3152 
3153         float tstFraction = 0.f;
3154         if (!hitWithNoTS) {
3155           const auto findTSIt =
3156               std::find(detIdToTracksterId_Map[hitDetId].begin(),
3157                         detIdToTracksterId_Map[hitDetId].end(),
3158                         HGVHistoProducerAlgo::detIdInfoInTrackster{
3159                             tstId, 0, 0.f});  // only the first element is used for the matching (overloaded operator==)
3160           if (findTSIt != detIdToTracksterId_Map[hitDetId].end()) {
3161             if (valType == 0) {
3162               tstFraction = findTSIt->fraction;
3163             } else {
3164               const auto iLC = std::find(
3165                   tracksters[tstId].vertices().begin(), tracksters[tstId].vertices().end(), findTSIt->clusterId);
3166               if (iLC != tracksters[tstId].vertices().end()) {
3167                 tstFraction = 1.f / tracksters[tstId].vertex_multiplicity(
3168                                         std::distance(std::begin(tracksters[tstId].vertices()), iLC));
3169               }
3170             }
3171           }
3172         }
3173         // Here do not divide as before by the trackster energy weight. Should sum first
3174         // over all layers and divide with the total CP energy over all layers.
3175         if (tsPair.second.second == FLT_MAX) {
3176           tsPair.second.second = 0.f;
3177         }
3178         tsPair.second.second += min(pow(tstFraction - cpFraction, 2), pow(cpFraction, 2)) * hitEnergyWeight;
3179 
3180         LogDebug("HGCalValidator") << "\nTracksterId:\t" << tstId << "\tSimTracksterId:\t" << iSTS << "\tcpId:\t"
3181                                    << cpId << "\ttstfraction, cpfraction:\t" << tstFraction << ", " << cpFraction
3182                                    << "\thitEnergyWeight:\t" << hitEnergyWeight << "\tadded delta:\t"
3183                                    << pow((tstFraction - cpFraction), 2) * hitEnergyWeight
3184                                    << "\tcurrent Sim-score numerator:\t" << tsPair.second.second
3185                                    << "\tshared Sim energy:\t" << tsPair.second.first << '\n';
3186       }
3187     }  // end of loop through SimCluster SimHits on current layer
3188 
3189     if (simOnLayer.layerClusterIdToEnergyAndScore.empty())
3190       LogDebug("HGCalValidator") << "CP Id:\t" << cpId << "\tTS id:\t-1"
3191                                  << " Sub score in \t -1\n";
3192 
3193     for (const auto& tsPair : simOnLayer.layerClusterIdToEnergyAndScore) {
3194       const auto tstId = tsPair.first;
3195       // 3D score here without the denominator at this point
3196       if (score3d_iSTS[tstId] == FLT_MAX) {
3197         score3d_iSTS[tstId] = 0.f;
3198       }
3199       score3d_iSTS[tstId] += tsPair.second.second;
3200       tstSharedEnergy[iSTS][tstId] += tsPair.second.first;
3201     }
3202     //} // end of loop through layers
3203 
3204     const auto scoreDenom = (valType == 0) ? SimEnergyWeight : hitsEnergyWeight;
3205     const auto energyDenom = (valType == 0) ? SimEnergy : SimEnergy_LC;
3206 
3207     const auto sts_eta = sts.barycenter().eta();
3208     const auto sts_phi = sts.barycenter().phi();
3209     const auto sts_en = sts.raw_energy();
3210     const auto sts_pt = sts.raw_pt();
3211     histograms.h_denom_caloparticle_eta[valType][count]->Fill(sts_eta);
3212     histograms.h_denom_caloparticle_phi[valType][count]->Fill(sts_phi);
3213     histograms.h_denom_caloparticle_en[valType][count]->Fill(sts_en);
3214     histograms.h_denom_caloparticle_pt[valType][count]->Fill(sts_pt);
3215 
3216     //Loop through related Tracksters here
3217     // In case the threshold to associate a CaloParticle to a Trackster is
3218     // below 50%, there could be cases in which the CP is linked to more than
3219     // one tracksters, leading to efficiencies >1. This boolean is used to
3220     // avoid "over counting".
3221     bool sts_considered_efficient = false;
3222     bool sts_considered_pure = false;
3223     for (const auto tstId : stsId_tstId_related) {
3224       // Now time for the denominator
3225       score3d_iSTS[tstId] /= scoreDenom;
3226       const auto tstSharedEnergyFrac = tstSharedEnergy[iSTS][tstId] / energyDenom;
3227       LogDebug("HGCalValidator") << "STS id: " << iSTS << "\t(CP id: " << cpId << ")\tTS id: " << tstId
3228                                  << "\nSimEnergy: " << energyDenom << "\tSimEnergyWeight: " << SimEnergyWeight
3229                                  << "\tTrackste energy: " << tracksters[tstId].raw_energy()
3230                                  << "\nscore: " << score3d_iSTS[tstId]
3231                                  << "\tshared energy: " << tstSharedEnergy[iSTS][tstId]
3232                                  << "\tshared energy fraction: " << tstSharedEnergyFrac << "\n";
3233 
3234       histograms.h_score_caloparticle2trackster[valType][count]->Fill(score3d_iSTS[tstId]);
3235       histograms.h_sharedenergy_caloparticle2trackster[valType][count]->Fill(tstSharedEnergyFrac);
3236       histograms.h_energy_vs_score_caloparticle2trackster[valType][count]->Fill(score3d_iSTS[tstId],
3237                                                                                 tstSharedEnergyFrac);
3238       // Fill the numerator for the efficiency calculation. The efficiency is computed by considering the energy shared between a Trackster and a _corresponding_ caloParticle. The threshold is configurable via python.
3239       if (!sts_considered_efficient && (tstSharedEnergyFrac >= minTSTSharedEneFracEfficiency_)) {
3240         sts_considered_efficient = true;
3241         histograms.h_numEff_caloparticle_eta[valType][count]->Fill(sts_eta);
3242         histograms.h_numEff_caloparticle_phi[valType][count]->Fill(sts_phi);
3243         histograms.h_numEff_caloparticle_en[valType][count]->Fill(sts_en);
3244         histograms.h_numEff_caloparticle_pt[valType][count]->Fill(sts_pt);
3245       }
3246 
3247       if (score3d_iSTS[tstId] < ScoreCutSTStoTSPurDup) {
3248         if (tracksters_PurityDuplicate[tstId] < 1)
3249           tracksters_PurityDuplicate[tstId]++;  // for Purity
3250         if (sts_considered_pure)
3251           tracksters_PurityDuplicate[tstId]++;  // for Duplicate
3252         sts_considered_pure = true;
3253       }
3254     }  // end of loop through Tracksters related to SimTrackster
3255 
3256     const auto best = std::min_element(std::begin(score3d_iSTS), std::end(score3d_iSTS));
3257     if (best != score3d_iSTS.end()) {
3258       const auto bestTstId = std::distance(std::begin(score3d_iSTS), best);
3259       const auto bestTstSharedEnergyFrac = tstSharedEnergy[iSTS][bestTstId] / energyDenom;
3260       histograms.h_scorePur_caloparticle2trackster[valType][count]->Fill(*best);
3261       histograms.h_sharedenergy_caloparticle2trackster_assoc[valType][count]->Fill(bestTstSharedEnergyFrac);
3262       histograms.h_sharedenergy_caloparticle2trackster_assoc_vs_eta[valType][count]->Fill(sts_eta,
3263                                                                                           bestTstSharedEnergyFrac);
3264       histograms.h_sharedenergy_caloparticle2trackster_assoc_vs_phi[valType][count]->Fill(sts_phi,
3265                                                                                           bestTstSharedEnergyFrac);
3266       histograms.h_energy_vs_score_caloparticle2bestTrackster[valType][count]->Fill(*best, bestTstSharedEnergyFrac);
3267       LogDebug("HGCalValidator") << count << " " << sts_eta << " " << sts_phi << " "
3268                                  << tracksters[bestTstId].raw_energy() << " " << sts.raw_energy() << " "
3269                                  << bestTstSharedEnergyFrac << "\n";
3270 
3271       if (score3d_iSTS.size() > 1) {
3272         auto best2 = (best == score3d_iSTS.begin()) ? std::next(best, 1) : score3d_iSTS.begin();
3273         for (auto tstId = score3d_iSTS.begin(); tstId != score3d_iSTS.end() && tstId != best; tstId++)
3274           if (*tstId < *best2)
3275             best2 = tstId;
3276         const auto best2TstId = std::distance(std::begin(score3d_iSTS), best2);
3277         const auto best2TstSharedEnergyFrac = tstSharedEnergy[iSTS][best2TstId] / energyDenom;
3278         histograms.h_scoreDupl_caloparticle2trackster[valType][count]->Fill(*best2);
3279         histograms.h_sharedenergy_caloparticle2trackster_assoc2[valType][count]->Fill(best2TstSharedEnergyFrac);
3280         histograms.h_energy_vs_score_caloparticle2bestTrackster2[valType][count]->Fill(*best2,
3281                                                                                        best2TstSharedEnergyFrac);
3282       }
3283     }
3284   }  // end of loop through SimTracksters
3285 
3286   // Fill the plots to compute the different metrics linked to
3287   // reco-level, namely fake-rate an merge-rate. Should *not*
3288   // restrict only to the selected caloParaticles.
3289   for (unsigned int tstId = 0; tstId < nTracksters; ++tstId) {
3290     const auto& tst = tracksters[tstId];
3291     if (tst.vertices().empty())
3292       continue;
3293     const auto iTS_eta = tst.barycenter().eta();
3294     const auto iTS_phi = tst.barycenter().phi();
3295     const auto iTS_en = tst.raw_energy();
3296     const auto iTS_pt = tst.raw_pt();
3297     histograms.h_denom_trackster_eta[valType][count]->Fill(iTS_eta);
3298     histograms.h_denom_trackster_phi[valType][count]->Fill(iTS_phi);
3299     histograms.h_denom_trackster_en[valType][count]->Fill(iTS_en);
3300     histograms.h_denom_trackster_pt[valType][count]->Fill(iTS_pt);
3301 
3302     if (tracksters_PurityDuplicate[tstId] > 0) {
3303       histograms.h_num_caloparticle_eta[valType][count]->Fill(iTS_eta);
3304       histograms.h_num_caloparticle_phi[valType][count]->Fill(iTS_phi);
3305       histograms.h_num_caloparticle_en[valType][count]->Fill(iTS_en);
3306       histograms.h_num_caloparticle_pt[valType][count]->Fill(iTS_pt);
3307 
3308       if (tracksters_PurityDuplicate[tstId] > 1) {
3309         histograms.h_numDup_trackster_eta[valType][count]->Fill(iTS_eta);
3310         histograms.h_numDup_trackster_phi[valType][count]->Fill(iTS_phi);
3311         histograms.h_numDup_trackster_en[valType][count]->Fill(iTS_en);
3312         histograms.h_numDup_trackster_pt[valType][count]->Fill(iTS_pt);
3313       }
3314     }
3315 
3316     if (tracksters_FakeMerge[tstId] > 0) {
3317       histograms.h_num_trackster_eta[valType][count]->Fill(iTS_eta);
3318       histograms.h_num_trackster_phi[valType][count]->Fill(iTS_phi);
3319       histograms.h_num_trackster_en[valType][count]->Fill(iTS_en);
3320       histograms.h_num_trackster_pt[valType][count]->Fill(iTS_pt);
3321 
3322       if (tracksters_FakeMerge[tstId] > 1) {
3323         histograms.h_numMerge_trackster_eta[valType][count]->Fill(iTS_eta);
3324         histograms.h_numMerge_trackster_phi[valType][count]->Fill(iTS_phi);
3325         histograms.h_numMerge_trackster_en[valType][count]->Fill(iTS_en);
3326         histograms.h_numMerge_trackster_pt[valType][count]->Fill(iTS_pt);
3327       }
3328     }
3329   }  // End loop over Tracksters
3330 }
3331 
3332 void HGVHistoProducerAlgo::fill_trackster_histos(
3333     const Histograms& histograms,
3334     const int count,
3335     const ticl::TracksterCollection& tracksters,
3336     const reco::CaloClusterCollection& layerClusters,
3337     const ticl::TracksterCollection& simTSs,
3338     const ticl::TracksterCollection& simTSs_fromCP,
3339     const std::map<unsigned int, std::vector<unsigned int>>& cpToSc_SimTrackstersMap,
3340     std::vector<SimCluster> const& sC,
3341     const edm::ProductID& cPHandle_id,
3342     std::vector<CaloParticle> const& cP,
3343     std::vector<size_t> const& cPIndices,
3344     std::vector<size_t> const& cPSelectedIndices,
3345     std::unordered_map<DetId, const HGCRecHit*> const& hitMap,
3346     unsigned int layers) const {
3347   //Each event to be treated as two events:
3348   //an event in +ve endcap, plus another event in -ve endcap.
3349 
3350   //To keep track of total num of Tracksters
3351   int totNTstZm = 0;  //-z
3352   int totNTstZp = 0;  //+z
3353   //To count the number of Tracksters with 3 contiguous layers per event.
3354   int totNContTstZp = 0;  //+z
3355   int totNContTstZm = 0;  //-z
3356   //For the number of Tracksters without 3 contiguous layers per event.
3357   int totNNotContTstZp = 0;  //+z
3358   int totNNotContTstZm = 0;  //-z
3359   // Check below the score of cont and non cont Tracksters
3360   std::vector<bool> contTracksters;
3361   contTracksters.clear();
3362 
3363   //[tstId]-> vector of 2d layer clusters size
3364   std::unordered_map<unsigned int, std::vector<unsigned int>> multiplicity;
3365   //[tstId]-> [layer][cluster size]
3366   std::unordered_map<unsigned int, std::vector<unsigned int>> multiplicity_vs_layer;
3367   //We will need for the scale text option
3368   // unsigned int totalLcInTsts = 0;
3369   // for (unsigned int tstId = 0; tstId < nTracksters; ++tstId) {
3370   //   totalLcInTsts = totalLcInTsts + tracksters[tstId].vertices().size();
3371   // }
3372 
3373   const auto nTracksters = tracksters.size();
3374   // loop through Tracksters
3375   for (unsigned int tstId = 0; tstId < nTracksters; ++tstId) {
3376     const auto& tst = tracksters[tstId];
3377     if (tst.vertices().empty())
3378       continue;
3379 
3380     if (tst.barycenter().z() < 0.)
3381       totNTstZm++;
3382     else if (tst.barycenter().z() > 0.)
3383       totNTstZp++;
3384 
3385     //Total number of layer clusters in Trackster
3386     int tnLcInTst = 0;
3387 
3388     //To keep track of total num of layer clusters per Trackster
3389     //tnLcInTstperlaypz[layerid], tnLcInTstperlaymz[layerid]
3390     std::vector<int> tnLcInTstperlay(1000, 0);  //+z
3391 
3392     //For the layers the Trackster expands to. Will use a set because there would be many
3393     //duplicates and then go back to vector for random access, since they say it is faster.
3394     std::set<unsigned int> trackster_layers;
3395 
3396     bool tracksterInZplus = false;
3397     bool tracksterInZminus = false;
3398 
3399     //Loop through layer clusters
3400     for (const auto lcId : tst.vertices()) {
3401       //take the hits and their fraction of the specific layer cluster.
3402       const auto& hits_and_fractions = layerClusters[lcId].hitsAndFractions();
3403 
3404       //For the multiplicity of the 2d layer clusters in Tracksters
3405       multiplicity[tstId].emplace_back(hits_and_fractions.size());
3406 
3407       const auto firstHitDetId = hits_and_fractions[0].first;
3408       //The layer that the layer cluster belongs to
3409       const auto layerid = recHitTools_->getLayerWithOffset(firstHitDetId) +
3410                            layers * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
3411       trackster_layers.insert(layerid);
3412       multiplicity_vs_layer[tstId].emplace_back(layerid);
3413 
3414       tnLcInTstperlay[layerid]++;
3415       tnLcInTst++;
3416 
3417       if (recHitTools_->zside(firstHitDetId) > 0.)
3418         tracksterInZplus = true;
3419       else if (recHitTools_->zside(firstHitDetId) < 0.)
3420         tracksterInZminus = true;
3421     }  // end of loop through layerClusters
3422 
3423     // Per layer : Loop 0->99
3424     for (unsigned ilayer = 0; ilayer < layers * 2; ++ilayer) {
3425       if (histograms.h_clusternum_in_trackster_perlayer[count].count(ilayer) && tnLcInTstperlay[ilayer] != 0) {
3426         histograms.h_clusternum_in_trackster_perlayer[count].at(ilayer)->Fill((float)tnLcInTstperlay[ilayer]);
3427       }
3428       // For the profile now of 2d layer cluster in Tracksters vs layer number.
3429       if (tnLcInTstperlay[ilayer] != 0) {
3430         histograms.h_clusternum_in_trackster_vs_layer[count]->Fill((float)ilayer, (float)tnLcInTstperlay[ilayer]);
3431       }
3432     }  // end of loop over layers
3433 
3434     // Looking for Tracksters with 3 contiguous layers per event.
3435     std::vector<int> trackster_layers_vec(trackster_layers.begin(), trackster_layers.end());
3436     // Check also for non contiguous Tracksters
3437     bool contiTrackster = false;
3438     // Start from 1 and go up to size - 1 element.
3439     if (trackster_layers_vec.size() >= 3) {
3440       for (unsigned int iLayer = 1; iLayer < trackster_layers_vec.size() - 1; ++iLayer) {
3441         if ((trackster_layers_vec[iLayer - 1] + 1 == trackster_layers_vec[iLayer]) &&
3442             (trackster_layers_vec[iLayer + 1] - 1 == trackster_layers_vec[iLayer])) {
3443           // Trackster with 3 contiguous layers per event
3444           if (tracksterInZplus)
3445             totNContTstZp++;
3446           else if (tracksterInZminus)
3447             totNContTstZm++;
3448 
3449           contiTrackster = true;
3450           break;
3451         }
3452       }
3453     }
3454     // Count non contiguous Tracksters
3455     if (!contiTrackster) {
3456       if (tracksterInZplus)
3457         totNNotContTstZp++;
3458       else if (tracksterInZminus)
3459         totNNotContTstZm++;
3460     }
3461 
3462     // Save for the score
3463     contTracksters.push_back(contiTrackster);
3464 
3465     histograms.h_clusternum_in_trackster[count]->Fill(tnLcInTst);
3466 
3467     for (unsigned int lc = 0; lc < multiplicity[tstId].size(); ++lc) {
3468       //multiplicity of the current LC
3469       float mlp = std::count(std::begin(multiplicity[tstId]), std::end(multiplicity[tstId]), multiplicity[tstId][lc]);
3470       //LogDebug("HGCalValidator") << "mlp %" << (100. * mlp)/ ((float) nLayerClusters) << std::endl;
3471       // histograms.h_multiplicityOfLCinTST[count]->Fill( mlp , multiplicity[tstId][lc] , 100. / (float) totalLcInTsts );
3472       histograms.h_multiplicityOfLCinTST[count]->Fill(mlp, multiplicity[tstId][lc]);
3473       //When plotting with the text option we want the entries to be the same
3474       //as the % of the current cell over the whole number of layerClusters. For this we need an extra histo.
3475       histograms.h_multiplicity_numberOfEventsHistogram[count]->Fill(mlp);
3476       //For the cluster multiplicity vs layer
3477       //First with the -z endcap (V10:0->49)
3478       if (multiplicity_vs_layer[tstId][lc] < layers) {
3479         histograms.h_multiplicityOfLCinTST_vs_layercluster_zminus[count]->Fill(mlp, multiplicity_vs_layer[tstId][lc]);
3480         histograms.h_multiplicity_zminus_numberOfEventsHistogram[count]->Fill(mlp);
3481       } else {  //Then for the +z (V10:50->99)
3482         histograms.h_multiplicityOfLCinTST_vs_layercluster_zplus[count]->Fill(
3483             mlp, multiplicity_vs_layer[tstId][lc] - layers);
3484         histograms.h_multiplicity_zplus_numberOfEventsHistogram[count]->Fill(mlp);
3485       }
3486       //For the cluster multiplicity vs cluster energy
3487       histograms.h_multiplicityOfLCinTST_vs_layerclusterenergy[count]->Fill(mlp,
3488                                                                             layerClusters[tst.vertices(lc)].energy());
3489     }
3490 
3491     if (!trackster_layers.empty()) {
3492       histograms.h_trackster_x[count]->Fill(tst.barycenter().x());
3493       histograms.h_trackster_y[count]->Fill(tst.barycenter().y());
3494       histograms.h_trackster_z[count]->Fill(tst.barycenter().z());
3495       histograms.h_trackster_eta[count]->Fill(tst.barycenter().eta());
3496       histograms.h_trackster_phi[count]->Fill(tst.barycenter().phi());
3497 
3498       histograms.h_trackster_firstlayer[count]->Fill((float)*trackster_layers.begin());
3499       histograms.h_trackster_lastlayer[count]->Fill((float)*trackster_layers.rbegin());
3500       histograms.h_trackster_layersnum[count]->Fill((float)trackster_layers.size());
3501 
3502       histograms.h_trackster_pt[count]->Fill(tst.raw_pt());
3503       histograms.h_trackster_energy[count]->Fill(tst.raw_energy());
3504     }
3505 
3506   }  //end of loop through Tracksters
3507 
3508   histograms.h_tracksternum[count]->Fill(totNTstZm + totNTstZp);
3509   histograms.h_conttracksternum[count]->Fill(totNContTstZp + totNContTstZm);
3510   histograms.h_nonconttracksternum[count]->Fill(totNNotContTstZp + totNNotContTstZm);
3511 
3512   // CaloParticle
3513   tracksters_to_SimTracksters(histograms,
3514                               count,
3515                               tracksters,
3516                               layerClusters,
3517                               simTSs_fromCP,
3518                               Linking,
3519                               simTSs_fromCP,
3520                               cpToSc_SimTrackstersMap,
3521                               sC,
3522                               cPHandle_id,
3523                               cP,
3524                               cPIndices,
3525                               cPSelectedIndices,
3526                               hitMap,
3527                               layers);
3528 
3529   // Pattern recognition
3530   tracksters_to_SimTracksters(histograms,
3531                               count,
3532                               tracksters,
3533                               layerClusters,
3534                               simTSs,
3535                               PatternRecognition,
3536                               simTSs_fromCP,
3537                               cpToSc_SimTrackstersMap,
3538                               sC,
3539                               cPHandle_id,
3540                               cP,
3541                               cPIndices,
3542                               cPSelectedIndices,
3543                               hitMap,
3544                               layers);
3545 }
3546 
3547 double HGVHistoProducerAlgo::distance2(const double x1,
3548                                        const double y1,
3549                                        const double x2,
3550                                        const double y2) const {  //distance squared
3551   const double dx = x1 - x2;
3552   const double dy = y1 - y2;
3553   return (dx * dx + dy * dy);
3554 }  //distance squaredq
3555 double HGVHistoProducerAlgo::distance(const double x1,
3556                                       const double y1,
3557                                       const double x2,
3558                                       const double y2) const {  //2-d distance on the layer (x-y)
3559   return std::sqrt(distance2(x1, y1, x2, y2));
3560 }
3561 
3562 void HGVHistoProducerAlgo::setRecHitTools(std::shared_ptr<hgcal::RecHitTools> recHitTools) {
3563   recHitTools_ = recHitTools;
3564 }
3565 
3566 DetId HGVHistoProducerAlgo::findmaxhit(const reco::CaloCluster& cluster,
3567                                        std::unordered_map<DetId, const HGCRecHit*> const& hitMap) const {
3568   const auto& hits_and_fractions = cluster.hitsAndFractions();
3569 
3570   DetId themaxid;
3571   double maxene = 0.;
3572   for (std::vector<std::pair<DetId, float>>::const_iterator it_haf = hits_and_fractions.begin();
3573        it_haf != hits_and_fractions.end();
3574        ++it_haf) {
3575     const DetId rh_detid = it_haf->first;
3576     const auto hitEn = hitMap.find(rh_detid)->second->energy();
3577 
3578     if (maxene < hitEn) {
3579       maxene = hitEn;
3580       themaxid = rh_detid;
3581     }
3582   }
3583 
3584   return themaxid;
3585 }
3586 
3587 double HGVHistoProducerAlgo::getEta(double eta) const {
3588   if (useFabsEta_)
3589     return fabs(eta);
3590   else
3591     return eta;
3592 }