Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-03-30 23:47:47

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