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
0014
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};
0020 const double ScoreCutSTStoTSPurDup_[] = {0.2, FLT_MIN};
0021
0022 HGVHistoProducerAlgo::HGVHistoProducerAlgo(const edm::ParameterSet& pset)
0023 :
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
0030 minEne_(pset.getParameter<double>("minEne")),
0031 maxEne_(pset.getParameter<double>("maxEne")),
0032 nintEne_(pset.getParameter<int>("nintEne")),
0033
0034
0035 minPt_(pset.getParameter<double>("minPt")),
0036 maxPt_(pset.getParameter<double>("maxPt")),
0037 nintPt_(pset.getParameter<int>("nintPt")),
0038
0039
0040 minPhi_(pset.getParameter<double>("minPhi")),
0041 maxPhi_(pset.getParameter<double>("maxPhi")),
0042 nintPhi_(pset.getParameter<int>("nintPhi")),
0043
0044
0045 minMixedHitsSimCluster_(pset.getParameter<double>("minMixedHitsSimCluster")),
0046 maxMixedHitsSimCluster_(pset.getParameter<double>("maxMixedHitsSimCluster")),
0047 nintMixedHitsSimCluster_(pset.getParameter<int>("nintMixedHitsSimCluster")),
0048
0049
0050 minMixedHitsCluster_(pset.getParameter<double>("minMixedHitsCluster")),
0051 maxMixedHitsCluster_(pset.getParameter<double>("maxMixedHitsCluster")),
0052 nintMixedHitsCluster_(pset.getParameter<int>("nintMixedHitsCluster")),
0053
0054
0055 minEneCl_(pset.getParameter<double>("minEneCl")),
0056 maxEneCl_(pset.getParameter<double>("maxEneCl")),
0057 nintEneCl_(pset.getParameter<int>("nintEneCl")),
0058
0059
0060 minLongDepBary_(pset.getParameter<double>("minLongDepBary")),
0061 maxLongDepBary_(pset.getParameter<double>("maxLongDepBary")),
0062 nintLongDepBary_(pset.getParameter<int>("nintLongDepBary")),
0063
0064
0065 minZpos_(pset.getParameter<double>("minZpos")),
0066 maxZpos_(pset.getParameter<double>("maxZpos")),
0067 nintZpos_(pset.getParameter<int>("nintZpos")),
0068
0069
0070 minTotNsimClsperlay_(pset.getParameter<double>("minTotNsimClsperlay")),
0071 maxTotNsimClsperlay_(pset.getParameter<double>("maxTotNsimClsperlay")),
0072 nintTotNsimClsperlay_(pset.getParameter<int>("nintTotNsimClsperlay")),
0073
0074
0075 minTotNClsperlay_(pset.getParameter<double>("minTotNClsperlay")),
0076 maxTotNClsperlay_(pset.getParameter<double>("maxTotNClsperlay")),
0077 nintTotNClsperlay_(pset.getParameter<int>("nintTotNClsperlay")),
0078
0079
0080 minEneClperlay_(pset.getParameter<double>("minEneClperlay")),
0081 maxEneClperlay_(pset.getParameter<double>("maxEneClperlay")),
0082 nintEneClperlay_(pset.getParameter<int>("nintEneClperlay")),
0083
0084
0085
0086
0087 minScore_(pset.getParameter<double>("minScore")),
0088 maxScore_(pset.getParameter<double>("maxScore")),
0089 nintScore_(pset.getParameter<int>("nintScore")),
0090
0091
0092
0093
0094
0095
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
0102 minTSTSharedEneFrac_(pset.getParameter<double>("minTSTSharedEneFrac")),
0103 maxTSTSharedEneFrac_(pset.getParameter<double>("maxTSTSharedEneFrac")),
0104 nintTSTSharedEneFrac_(pset.getParameter<int>("nintTSTSharedEneFrac")),
0105
0106
0107 minTotNsimClsperthick_(pset.getParameter<double>("minTotNsimClsperthick")),
0108 maxTotNsimClsperthick_(pset.getParameter<double>("maxTotNsimClsperthick")),
0109 nintTotNsimClsperthick_(pset.getParameter<int>("nintTotNsimClsperthick")),
0110
0111
0112 minTotNClsperthick_(pset.getParameter<double>("minTotNClsperthick")),
0113 maxTotNClsperthick_(pset.getParameter<double>("maxTotNClsperthick")),
0114 nintTotNClsperthick_(pset.getParameter<int>("nintTotNClsperthick")),
0115
0116
0117 minTotNcellsperthickperlayer_(pset.getParameter<double>("minTotNcellsperthickperlayer")),
0118 maxTotNcellsperthickperlayer_(pset.getParameter<double>("maxTotNcellsperthickperlayer")),
0119 nintTotNcellsperthickperlayer_(pset.getParameter<int>("nintTotNcellsperthickperlayer")),
0120
0121
0122 minDisToSeedperthickperlayer_(pset.getParameter<double>("minDisToSeedperthickperlayer")),
0123 maxDisToSeedperthickperlayer_(pset.getParameter<double>("maxDisToSeedperthickperlayer")),
0124 nintDisToSeedperthickperlayer_(pset.getParameter<int>("nintDisToSeedperthickperlayer")),
0125
0126
0127 minDisToSeedperthickperlayerenewei_(pset.getParameter<double>("minDisToSeedperthickperlayerenewei")),
0128 maxDisToSeedperthickperlayerenewei_(pset.getParameter<double>("maxDisToSeedperthickperlayerenewei")),
0129 nintDisToSeedperthickperlayerenewei_(pset.getParameter<int>("nintDisToSeedperthickperlayerenewei")),
0130
0131
0132 minDisToMaxperthickperlayer_(pset.getParameter<double>("minDisToMaxperthickperlayer")),
0133 maxDisToMaxperthickperlayer_(pset.getParameter<double>("maxDisToMaxperthickperlayer")),
0134 nintDisToMaxperthickperlayer_(pset.getParameter<int>("nintDisToMaxperthickperlayer")),
0135
0136
0137 minDisToMaxperthickperlayerenewei_(pset.getParameter<double>("minDisToMaxperthickperlayerenewei")),
0138 maxDisToMaxperthickperlayerenewei_(pset.getParameter<double>("maxDisToMaxperthickperlayerenewei")),
0139 nintDisToMaxperthickperlayerenewei_(pset.getParameter<int>("nintDisToMaxperthickperlayerenewei")),
0140
0141
0142 minDisSeedToMaxperthickperlayer_(pset.getParameter<double>("minDisSeedToMaxperthickperlayer")),
0143 maxDisSeedToMaxperthickperlayer_(pset.getParameter<double>("maxDisSeedToMaxperthickperlayer")),
0144 nintDisSeedToMaxperthickperlayer_(pset.getParameter<int>("nintDisSeedToMaxperthickperlayer")),
0145
0146
0147 minClEneperthickperlayer_(pset.getParameter<double>("minClEneperthickperlayer")),
0148 maxClEneperthickperlayer_(pset.getParameter<double>("maxClEneperthickperlayer")),
0149 nintClEneperthickperlayer_(pset.getParameter<int>("nintClEneperthickperlayer")),
0150
0151
0152 minCellsEneDensperthick_(pset.getParameter<double>("minCellsEneDensperthick")),
0153 maxCellsEneDensperthick_(pset.getParameter<double>("maxCellsEneDensperthick")),
0154 nintCellsEneDensperthick_(pset.getParameter<int>("nintCellsEneDensperthick")),
0155
0156
0157
0158 minTotNTSTs_(pset.getParameter<double>("minTotNTSTs")),
0159 maxTotNTSTs_(pset.getParameter<double>("maxTotNTSTs")),
0160 nintTotNTSTs_(pset.getParameter<int>("nintTotNTSTs")),
0161
0162
0163 minTotNClsinTSTs_(pset.getParameter<double>("minTotNClsinTSTs")),
0164 maxTotNClsinTSTs_(pset.getParameter<double>("maxTotNClsinTSTs")),
0165 nintTotNClsinTSTs_(pset.getParameter<int>("nintTotNClsinTSTs")),
0166
0167
0168 minTotNClsinTSTsperlayer_(pset.getParameter<double>("minTotNClsinTSTsperlayer")),
0169 maxTotNClsinTSTsperlayer_(pset.getParameter<double>("maxTotNClsinTSTsperlayer")),
0170 nintTotNClsinTSTsperlayer_(pset.getParameter<int>("nintTotNClsinTSTsperlayer")),
0171
0172
0173 minMplofLCs_(pset.getParameter<double>("minMplofLCs")),
0174 maxMplofLCs_(pset.getParameter<double>("maxMplofLCs")),
0175 nintMplofLCs_(pset.getParameter<int>("nintMplofLCs")),
0176
0177
0178 minSizeCLsinTSTs_(pset.getParameter<double>("minSizeCLsinTSTs")),
0179 maxSizeCLsinTSTs_(pset.getParameter<double>("maxSizeCLsinTSTs")),
0180 nintSizeCLsinTSTs_(pset.getParameter<int>("nintSizeCLsinTSTs")),
0181
0182
0183 minClEnepermultiplicity_(pset.getParameter<double>("minClEnepermultiplicity")),
0184 maxClEnepermultiplicity_(pset.getParameter<double>("maxClEnepermultiplicity")),
0185 nintClEnepermultiplicity_(pset.getParameter<int>("nintClEnepermultiplicity")),
0186
0187
0188 minX_(pset.getParameter<double>("minX")),
0189 maxX_(pset.getParameter<double>("maxX")),
0190 nintX_(pset.getParameter<int>("nintX")),
0191
0192
0193 minY_(pset.getParameter<double>("minY")),
0194 maxY_(pset.getParameter<double>("maxY")),
0195 nintY_(pset.getParameter<int>("nintY")),
0196
0197
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
0301 std::string istr2 = "";
0302
0303 if (ilayer < layers) {
0304 istr2 = std::to_string(ilayer + 1) + " in z-";
0305 } else {
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 }
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 }
0324
0325
0326
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
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
0398 std::string istr2 = "";
0399
0400 if (ilayer < layers) {
0401 istr2 = std::to_string(ilayer + 1) + " in z-";
0402 } else {
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 }
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
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
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
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
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
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
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
0664 std::string istr2 = "";
0665
0666 if (ilayer < layers) {
0667 istr2 = std::to_string(ilayer + 1) + " in z-";
0668 } else {
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
0706 std::string istr2 = "";
0707
0708 if (ilayer < layers) {
0709 istr2 = std::to_string(ilayer + 1) + " in z-";
0710 } else {
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
0874 std::string istr2 = "";
0875
0876 if (ilayer < layers) {
0877 istr2 = std::to_string(ilayer + 1) + " in z-";
0878 } else {
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
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
0907 std::string istr3 = "";
0908
0909 if (ilayer < layers) {
0910 istr3 = std::to_string(ilayer + 1) + " in z- ";
0911 } else {
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
0981 std::string istr2 = "";
0982
0983 if (ilayer < layers) {
0984 istr2 = std::to_string(ilayer + 1) + " in z-";
0985 } else {
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
1119 const string val[] = {"_Link", "_PR"};
1120
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
1218
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
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
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
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
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
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
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
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
1428
1429
1430
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
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 }
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
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
1571
1572
1573
1574
1575
1576
1577
1578 std::vector<int> tnscpl(1000, 0);
1579
1580
1581 std::map<std::string, int> tnscpthplus;
1582 tnscpthplus.clear();
1583 std::map<std::string, int> tnscpthminus;
1584 tnscpthminus.clear();
1585
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
1591 tnscpthplus.insert(std::pair<std::string, int>("mixed", 0));
1592 tnscpthminus.insert(std::pair<std::string, int>("mixed", 0));
1593
1594
1595 for (const auto& sc : simClusters) {
1596
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
1606 double thickness = 0.;
1607
1608 std::vector<int> occurenceSCinlayer(1000, 0);
1609
1610
1611 for (const auto& hAndF : sc.hits_and_fractions()) {
1612 const DetId sh_detid = hAndF.first;
1613
1614
1615 int layerid =
1616 recHitTools_->getLayerWithOffset(sh_detid) + layers * ((recHitTools_->zside(sh_detid) + 1) >> 1) - 1;
1617
1618 int zside = recHitTools_->zside(sh_detid);
1619
1620
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 {
1652 LogDebug("HGCalValidator")
1653 << " You are running a geometry that contains thicknesses different than the normal ones. "
1654 << "\n";
1655 }
1656
1657 }
1658
1659
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
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
1674 tnscpthminus[std::to_string((int)thickness)]++;
1675 }
1676
1677 }
1678
1679
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
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
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
1710
1711
1712
1713
1714
1715
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
1755
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});
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
1787
1788
1789
1790
1791
1792
1793
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
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
1817
1818
1819
1820
1821
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
1834
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 }
1845
1846 }
1847
1848
1849
1850
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())
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())
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 }
1914
1915
1916
1917
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
2010
2011
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
2021
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
2033
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
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
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())
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
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
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())
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 }
2088
2089
2090
2091
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
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
2191
2192
2193
2194
2195
2196
2197
2198 std::vector<int> tnlcpl(1000, 0);
2199
2200
2201 std::map<std::string, int> tnlcpthplus;
2202 tnlcpthplus.clear();
2203 std::map<std::string, int> tnlcpthminus;
2204 tnlcpthminus.clear();
2205
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
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
2227
2228 std::vector<double> tecpl(1000, 0.0);
2229
2230 std::vector<double> ldbar(1000, 0.0);
2231
2232
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
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
2251 double maxx = recHitTools_->getPosition(maxid).x();
2252 double maxy = recHitTools_->getPosition(maxid).y();
2253
2254
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
2264 double thickness = 0.;
2265
2266 int layerid = 0;
2267
2268
2269 int lay = 0;
2270
2271 std::string istr = "";
2272
2273 bool cluslay = true;
2274
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
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
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 {
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
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
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 }
2373
2374
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
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
2389 tnlcpthminus[std::to_string((int)thickness)]++;
2390 }
2391
2392
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
2414 if (histograms.h_cellsnum_perthickperlayer.count(istr)) {
2415 histograms.h_cellsnum_perthickperlayer.at(istr)->Fill(hits_and_fractions.size());
2416 }
2417
2418
2419 double distancebetseedandmax = distance(seedx, seedy, maxx, maxy);
2420
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
2433 tecpl[layerid] = tecpl[layerid] + lc_en;
2434 ldbar[layerid] = ldbar[layerid] + lc_en * cummatbudg[(double)lay];
2435
2436 }
2437
2438
2439 double sumeneallcluspl = 0.;
2440 double sumeneallclusmi = 0.;
2441
2442 double sumldbarpl = 0.;
2443 double sumldbarmi = 0.;
2444
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
2450
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
2458 sumeneallclusmi = sumeneallclusmi + tecpl[ilayer];
2459
2460 sumldbarmi = sumldbarmi + ldbar[ilayer];
2461 } else {
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
2468 sumeneallcluspl = sumeneallcluspl + tecpl[ilayer];
2469
2470 sumldbarpl = sumldbarpl + ldbar[ilayer];
2471 }
2472
2473 }
2474
2475
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
2483 histograms.h_mixedhitscluster_zplus[count]->Fill(tnlcpthplus["mixed"]);
2484 histograms.h_mixedhitscluster_zminus[count]->Fill(tnlcpthminus["mixed"]);
2485
2486
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
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
2524
2525
2526 std::vector<std::vector<std::pair<unsigned int, float>>> stsInTrackster;
2527 stsInTrackster.resize(nTracksters);
2528
2529
2530
2531
2532
2533 std::unordered_map<int, caloParticleOnLayer> cPOnLayer;
2534 std::unordered_map<int, std::vector<caloParticleOnLayer>> sCOnLayer;
2535
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 {
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())
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
2595 for (const auto& simCluster : cP[cpId].simClusters()) {
2596 auto iSim = simTSs[iSTS].seedIndex();
2597 if (simTSs[iSTS].seedID() != cPHandle_id) {
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
2607
2608 const auto itcheck = hitMap.find(hitid);
2609
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
2620
2621
2622
2623
2624
2625
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});
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
2644
2645
2646 cPOnLayer[cpId].energy += it_haf.second * hitEn;
2647 sCOnLayer[cpId][iSim].energy += elemFr * hitEn;
2648
2649
2650
2651
2652
2653
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
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 }
2672 }
2673 }
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);
2683 }
2684 });
2685 return hits_and_fractions_norm;
2686 };
2687
2688 auto ScoreCutSTStoTSPurDup = ScoreCutSTStoTSPurDup_[0];
2689 auto ScoreCutTStoSTSFakeMerge = ScoreCutTStoSTSFakeMerge_[0];
2690
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
2711
2712
2713
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
2723
2724
2725
2726
2727
2728
2729
2730
2731
2732
2733
2734
2735
2736
2737
2738 std::vector<int> hitsToCaloParticleId(numberOfHitsInTS);
2739
2740
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
2750
2751
2752
2753
2754
2755
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});
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
2776
2777
2778
2779
2780 if (rhFraction == 0.) {
2781 hitsToCaloParticleId[iHit] = -2;
2782 numberOfHaloHitsInTS++;
2783 }
2784
2785
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
2793 const auto hitEn = hitMap.find(rh_detid)->second->energy();
2794
2795
2796
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)
2806 iSim = 0;
2807
2808
2809
2810
2811
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
2818
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
2824
2825 stsInTrackster[tstId].emplace_back(iSTS, FLT_MAX);
2826
2827
2828 if (shared_fraction > maxCPEnergyInTS) {
2829
2830 maxCPEnergyInTS = CPEnergyInTS[cpId];
2831 maxCPId = cpId;
2832 }
2833 }
2834
2835 hitsToCaloParticleId[iHit] = maxCPId;
2836 }
2837
2838 }
2839
2840
2841
2842 for (auto c : hitsToCaloParticleId) {
2843 if (c < 0)
2844 numberOfNoiseHitsInTS++;
2845 else
2846 occurrencesCPinTS[c]++;
2847 }
2848
2849
2850
2851 for (auto& c : occurrencesCPinTS) {
2852 if (c.second > maxCPNumberOfHitsInTS) {
2853 maxCPId_byNumberOfHits = c.first;
2854 maxCPNumberOfHitsInTS = c.second;
2855 }
2856 }
2857
2858
2859 for (auto& c : CPEnergyInTS) {
2860 if (c.second > maxEnergySharedTSandCP) {
2861 maxCPId_byEnergy = c.first;
2862 maxEnergySharedTSandCP = c.second;
2863 }
2864 }
2865
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"
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 }
2891
2892
2893 for (unsigned int tstId = 0; tstId < nTracksters; ++tstId) {
2894 const auto& tst = tracksters[tstId];
2895 if (tst.vertices().empty())
2896 continue;
2897
2898
2899
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
2906 for (auto& stsPair : stsInTrackster[tstId]) {
2907
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
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});
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 }
2973
2974
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)
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 }
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 }
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
3047
3048
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)
3066 iSim = 0;
3067 auto& simOnLayer = (valType == 0) ? cPOnLayer[cpId] : sCOnLayer[cpId][iSim];
3068
3069
3070
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
3077 const auto SimNumberOfHits = simOnLayer.hits_and_fractions.size();
3078 if (SimNumberOfHits == 0)
3079 continue;
3080 SimEnergy += simOnLayer.energy;
3081 int tstWithMaxEnergyInCP = -1;
3082
3083 float maxEnergyTSperlayerinSim = 0.f;
3084 float SimEnergyFractionInTSperlayer = 0.f;
3085
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
3109
3110
3111
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;
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});
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
3157
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 }
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
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
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
3200
3201
3202
3203
3204 bool sts_considered_efficient = false;
3205 bool sts_considered_pure = false;
3206 for (const auto tstId : stsId_tstId_related) {
3207
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
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]++;
3233 if (sts_considered_pure)
3234 tracksters_PurityDuplicate[tstId]++;
3235 sts_considered_pure = true;
3236 }
3237 }
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 }
3268
3269
3270
3271
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 }
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
3331
3332
3333
3334 int totNTstZm = 0;
3335 int totNTstZp = 0;
3336
3337 int totNContTstZp = 0;
3338 int totNContTstZm = 0;
3339
3340 int totNNotContTstZp = 0;
3341 int totNNotContTstZm = 0;
3342
3343 std::vector<bool> contTracksters;
3344 contTracksters.clear();
3345
3346
3347 std::unordered_map<unsigned int, std::vector<unsigned int>> multiplicity;
3348
3349 std::unordered_map<unsigned int, std::vector<unsigned int>> multiplicity_vs_layer;
3350
3351
3352
3353
3354
3355
3356 const auto nTracksters = tracksters.size();
3357
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
3369 int tnLcInTst = 0;
3370
3371
3372
3373 std::vector<int> tnLcInTstperlay(1000, 0);
3374
3375
3376
3377 std::set<unsigned int> trackster_layers;
3378
3379 bool tracksterInZplus = false;
3380 bool tracksterInZminus = false;
3381
3382
3383 for (const auto lcId : tst.vertices()) {
3384
3385 const auto& hits_and_fractions = layerClusters[lcId].hitsAndFractions();
3386
3387
3388 multiplicity[tstId].emplace_back(hits_and_fractions.size());
3389
3390 const auto firstHitDetId = hits_and_fractions[0].first;
3391
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 }
3405
3406
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
3412 if (tnLcInTstperlay[ilayer] != 0) {
3413 histograms.h_clusternum_in_trackster_vs_layer[count]->Fill((float)ilayer, (float)tnLcInTstperlay[ilayer]);
3414 }
3415 }
3416
3417
3418 std::vector<int> trackster_layers_vec(trackster_layers.begin(), trackster_layers.end());
3419
3420 bool contiTrackster = false;
3421
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
3427 if (tracksterInZplus)
3428 totNContTstZp++;
3429 else if (tracksterInZminus)
3430 totNContTstZm++;
3431
3432 contiTrackster = true;
3433 break;
3434 }
3435 }
3436 }
3437
3438 if (!contiTrackster) {
3439 if (tracksterInZplus)
3440 totNNotContTstZp++;
3441 else if (tracksterInZminus)
3442 totNNotContTstZm++;
3443 }
3444
3445
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
3452 float mlp = std::count(std::begin(multiplicity[tstId]), std::end(multiplicity[tstId]), multiplicity[tstId][lc]);
3453
3454
3455 histograms.h_multiplicityOfLCinTST[count]->Fill(mlp, multiplicity[tstId][lc]);
3456
3457
3458 histograms.h_multiplicity_numberOfEventsHistogram[count]->Fill(mlp);
3459
3460
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 {
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
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 }
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
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
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 {
3534 const double dx = x1 - x2;
3535 const double dy = y1 - y2;
3536 return (dx * dx + dy * dy);
3537 }
3538 double HGVHistoProducerAlgo::distance(const double x1,
3539 const double y1,
3540 const double x2,
3541 const double y2) const {
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 }
