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