File indexing completed on 2025-07-09 05:00:40
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 rtos = ";score Reco-to-Sim";
1118 const string stor = ";score Sim-to-Reco";
1119 const string shREnFr = ";shared Reco energy fraction";
1120 const string shSEnFr = ";shared Sim energy fraction";
1121
1122 histograms.h_score_trackster2caloparticle[valType].push_back(
1123 ibook.book1D("Score_trackster2" + ref_[valType],
1124 "Score of Trackster per " + refText_[valType] + rtos,
1125 nintScore_,
1126 minScore_,
1127 maxScore_));
1128 histograms.h_score_trackster2bestCaloparticle[valType].push_back(
1129 ibook.book1D("ScoreFake_trackster2" + ref_[valType],
1130 "Score of Trackster per best " + refText_[valType] + rtos,
1131 nintScore_,
1132 minScore_,
1133 maxScore_));
1134 histograms.h_score_trackster2bestCaloparticle2[valType].push_back(
1135 ibook.book1D("ScoreMerge_trackster2" + ref_[valType],
1136 "Score of Trackster per 2^{nd} best " + refText_[valType] + rtos,
1137 nintScore_,
1138 minScore_,
1139 maxScore_));
1140 histograms.h_score_caloparticle2trackster[valType].push_back(
1141 ibook.book1D("Score_" + ref_[valType] + "2trackster",
1142 "Score of " + refText_[valType] + " per Trackster" + stor,
1143 nintScore_,
1144 minScore_,
1145 maxScore_));
1146 histograms.h_scorePur_caloparticle2trackster[valType].push_back(
1147 ibook.book1D("ScorePur_" + ref_[valType] + "2trackster",
1148 "Score of " + refText_[valType] + " per best Trackster" + stor,
1149 nintScore_,
1150 minScore_,
1151 maxScore_));
1152 histograms.h_scoreDupl_caloparticle2trackster[valType].push_back(
1153 ibook.book1D("ScoreDupl_" + ref_[valType] + "2trackster",
1154 "Score of " + refText_[valType] + " per 2^{nd} best Trackster" + stor,
1155 nintScore_,
1156 minScore_,
1157 maxScore_));
1158 histograms.h_energy_vs_score_trackster2caloparticle[valType].push_back(
1159 ibook.book2D("Energy_vs_Score_trackster2" + ref_[valType],
1160 "Energy vs Score of Trackster per " + refText_[valType] + rtos + shREnFr,
1161 nintScore_,
1162 minScore_,
1163 maxScore_,
1164 nintSharedEneFrac_,
1165 minTSTSharedEneFrac_,
1166 maxTSTSharedEneFrac_));
1167 histograms.h_energy_vs_score_trackster2bestCaloparticle[valType].push_back(
1168 ibook.book2D("Energy_vs_Score_trackster2best" + ref_[valType],
1169 "Energy vs Score of Trackster per best " + refText_[valType] + rtos + shREnFr,
1170 nintScore_,
1171 minScore_,
1172 maxScore_,
1173 nintSharedEneFrac_,
1174 minTSTSharedEneFrac_,
1175 maxTSTSharedEneFrac_));
1176 histograms.h_energy_vs_score_trackster2bestCaloparticle2[valType].push_back(
1177 ibook.book2D("Energy_vs_Score_trackster2secBest" + ref_[valType],
1178 "Energy vs Score of Trackster per 2^{nd} best " + refText_[valType] + rtos + shREnFr,
1179 nintScore_,
1180 minScore_,
1181 maxScore_,
1182 nintSharedEneFrac_,
1183 minTSTSharedEneFrac_,
1184 maxTSTSharedEneFrac_));
1185 histograms.h_energy_vs_score_caloparticle2trackster[valType].push_back(
1186 ibook.book2D("Energy_vs_Score_" + ref_[valType] + "2Trackster",
1187 "Energy vs Score of " + refText_[valType] + " per Trackster" + stor + shSEnFr,
1188 nintScore_,
1189 minScore_,
1190 maxScore_,
1191 nintSharedEneFrac_,
1192 minTSTSharedEneFrac_,
1193 maxTSTSharedEneFrac_));
1194 histograms.h_energy_vs_score_caloparticle2bestTrackster[valType].push_back(
1195 ibook.book2D("Energy_vs_Score_" + ref_[valType] + "2bestTrackster",
1196 "Energy vs Score of " + refText_[valType] + " per best Trackster" + stor + shSEnFr,
1197 nintScore_,
1198 minScore_,
1199 maxScore_,
1200 nintSharedEneFrac_,
1201 minTSTSharedEneFrac_,
1202 maxTSTSharedEneFrac_));
1203 histograms.h_energy_vs_score_caloparticle2bestTrackster2[valType].push_back(
1204 ibook.book2D("Energy_vs_Score_" + ref_[valType] + "2secBestTrackster",
1205 "Energy vs Score of " + refText_[valType] + " per 2^{nd} best Trackster" + stor + shSEnFr,
1206 nintScore_,
1207 minScore_,
1208 maxScore_,
1209 nintSharedEneFrac_,
1210 minTSTSharedEneFrac_,
1211 maxTSTSharedEneFrac_));
1212
1213
1214
1215 histograms.h_num_trackster_eta[valType].push_back(ibook.book1D(
1216 "Num_Trackster_Eta" + valSuffix_[valType], "Num Trackster Eta per Trackster;#eta", nintEta_, minEta_, maxEta_));
1217 histograms.h_numMerge_trackster_eta[valType].push_back(ibook.book1D("NumMerge_Trackster_Eta" + valSuffix_[valType],
1218 "Num Merge Trackster Eta per Trackster;#eta",
1219 nintEta_,
1220 minEta_,
1221 maxEta_));
1222 histograms.h_denom_trackster_eta[valType].push_back(ibook.book1D("Denom_Trackster_Eta" + valSuffix_[valType],
1223 "Denom Trackster Eta per Trackster;#eta",
1224 nintEta_,
1225 minEta_,
1226 maxEta_));
1227
1228 histograms.h_num_trackster_phi[valType].push_back(ibook.book1D(
1229 "Num_Trackster_Phi" + valSuffix_[valType], "Num Trackster Phi per Trackster;#phi", nintPhi_, minPhi_, maxPhi_));
1230 histograms.h_numMerge_trackster_phi[valType].push_back(ibook.book1D("NumMerge_Trackster_Phi" + valSuffix_[valType],
1231 "Num Merge Trackster Phi per Trackster;#phi",
1232 nintPhi_,
1233 minPhi_,
1234 maxPhi_));
1235 histograms.h_denom_trackster_phi[valType].push_back(ibook.book1D("Denom_Trackster_Phi" + valSuffix_[valType],
1236 "Denom Trackster Phi per Trackster;#phi",
1237 nintPhi_,
1238 minPhi_,
1239 maxPhi_));
1240
1241 histograms.h_num_trackster_en[valType].push_back(ibook.book1D("Num_Trackster_Energy" + valSuffix_[valType],
1242 "Num Trackster Energy per Trackster;energy [GeV]",
1243 nintEne_,
1244 minEne_,
1245 maxEne_));
1246 histograms.h_numMerge_trackster_en[valType].push_back(
1247 ibook.book1D("NumMerge_Trackster_Energy" + valSuffix_[valType],
1248 "Num Merge Trackster Energy per Trackster;energy [GeV]",
1249 nintEne_,
1250 minEne_,
1251 maxEne_));
1252 histograms.h_denom_trackster_en[valType].push_back(ibook.book1D("Denom_Trackster_Energy" + valSuffix_[valType],
1253 "Denom Trackster Energy per Trackster;energy [GeV]",
1254 nintEne_,
1255 minEne_,
1256 maxEne_));
1257
1258 histograms.h_num_trackster_pt[valType].push_back(ibook.book1D("Num_Trackster_Pt" + valSuffix_[valType],
1259 "Num Trackster p_{T} per Trackster;p_{T} [GeV]",
1260 nintPt_,
1261 minPt_,
1262 maxPt_));
1263 histograms.h_numMerge_trackster_pt[valType].push_back(
1264 ibook.book1D("NumMerge_Trackster_Pt" + valSuffix_[valType],
1265 "Num Merge Trackster p_{T} per Trackster;p_{T} [GeV]",
1266 nintPt_,
1267 minPt_,
1268 maxPt_));
1269 histograms.h_denom_trackster_pt[valType].push_back(ibook.book1D("Denom_Trackster_Pt" + valSuffix_[valType],
1270 "Denom Trackster p_{T} per Trackster;p_{T} [GeV]",
1271 nintPt_,
1272 minPt_,
1273 maxPt_));
1274
1275 histograms.h_sharedenergy_trackster2caloparticle[valType].push_back(
1276 ibook.book1D("SharedEnergy_trackster2" + ref_[valType],
1277 "Shared Energy of Trackster per " + refText_[valType] + shREnFr,
1278 nintSharedEneFrac_,
1279 minTSTSharedEneFrac_,
1280 maxTSTSharedEneFrac_));
1281 histograms.h_sharedenergy_trackster2bestCaloparticle[valType].push_back(
1282 ibook.book1D("SharedEnergy_trackster2" + ref_[valType] + "_assoc",
1283 "Shared Energy of Trackster per best " + refText_[valType] + shREnFr,
1284 nintSharedEneFrac_,
1285 minTSTSharedEneFrac_,
1286 maxTSTSharedEneFrac_));
1287 histograms.h_sharedenergy_trackster2bestCaloparticle_vs_eta[valType].push_back(ibook.bookProfile(
1288 "SharedEnergy_trackster2" + ref_[valType] + "_assoc_vs_eta",
1289 "Shared Energy of Trackster vs #eta per best " + refText_[valType] + ";Trackster #eta" + shREnFr,
1290 nintEta_,
1291 minEta_,
1292 maxEta_,
1293 minTSTSharedEneFrac_,
1294 maxTSTSharedEneFrac_));
1295 histograms.h_sharedenergy_trackster2bestCaloparticle_vs_phi[valType].push_back(ibook.bookProfile(
1296 "SharedEnergy_trackster2" + ref_[valType] + "_assoc_vs_phi",
1297 "Shared Energy of Trackster vs #phi per best " + refText_[valType] + ";Trackster #phi" + shREnFr,
1298 nintPhi_,
1299 minPhi_,
1300 maxPhi_,
1301 minTSTSharedEneFrac_,
1302 maxTSTSharedEneFrac_));
1303 histograms.h_sharedenergy_trackster2bestCaloparticle2[valType].push_back(
1304 ibook.book1D("SharedEnergy_trackster2" + ref_[valType] + "_assoc2",
1305 "Shared Energy of Trackster per 2^{nd} best " + refText_[valType] + shREnFr,
1306 nintSharedEneFrac_,
1307 minTSTSharedEneFrac_,
1308 maxTSTSharedEneFrac_));
1309
1310 histograms.h_sharedenergy_caloparticle2trackster[valType].push_back(
1311 ibook.book1D("SharedEnergy_" + ref_[valType] + "2trackster",
1312 "Shared Energy of " + refText_[valType] + " per Trackster" + shSEnFr,
1313 nintSharedEneFrac_,
1314 minTSTSharedEneFrac_,
1315 maxTSTSharedEneFrac_));
1316 histograms.h_sharedenergy_caloparticle2trackster_assoc[valType].push_back(
1317 ibook.book1D("SharedEnergy_" + ref_[valType] + "2trackster_assoc",
1318 "Shared Energy of " + refText_[valType] + " per best Trackster" + shSEnFr,
1319 nintSharedEneFrac_,
1320 minTSTSharedEneFrac_,
1321 maxTSTSharedEneFrac_));
1322 histograms.h_sharedenergy_caloparticle2trackster_assoc_vs_eta[valType].push_back(ibook.bookProfile(
1323 "SharedEnergy_" + ref_[valType] + "2trackster_assoc_vs_eta",
1324 "Shared Energy of " + refText_[valType] + " vs #eta per best Trackster;" + refText_[valType] + " #eta" + shSEnFr,
1325 nintEta_,
1326 minEta_,
1327 maxEta_,
1328 minTSTSharedEneFrac_,
1329 maxTSTSharedEneFrac_));
1330 histograms.h_sharedenergy_caloparticle2trackster_assoc_vs_phi[valType].push_back(ibook.bookProfile(
1331 "SharedEnergy_" + ref_[valType] + "2trackster_assoc_vs_phi",
1332 "Shared Energy of " + refText_[valType] + " vs #phi per best Trackster;" + refText_[valType] + " #phi" + shSEnFr,
1333 nintPhi_,
1334 minPhi_,
1335 maxPhi_,
1336 minTSTSharedEneFrac_,
1337 maxTSTSharedEneFrac_));
1338 histograms.h_sharedenergy_caloparticle2trackster_assoc2[valType].push_back(
1339 ibook.book1D("SharedEnergy_" + ref_[valType] + "2trackster_assoc2",
1340 "Shared Energy of " + refText_[valType] + " per 2^{nd} best Trackster;" + shSEnFr,
1341 nintSharedEneFrac_,
1342 minTSTSharedEneFrac_,
1343 maxTSTSharedEneFrac_));
1344
1345
1346 histograms.h_numEff_caloparticle_eta[valType].push_back(
1347 ibook.book1D("NumEff_" + ref_[valType] + "_Eta",
1348 "Num Efficiency " + refText_[valType] + " Eta per Trackster;#eta",
1349 nintEta_,
1350 minEta_,
1351 maxEta_));
1352 histograms.h_num_caloparticle_eta[valType].push_back(
1353 ibook.book1D("Num_" + ref_[valType] + "_Eta",
1354 "Num Purity " + refText_[valType] + " Eta per Trackster;#eta",
1355 nintEta_,
1356 minEta_,
1357 maxEta_));
1358 histograms.h_numDup_trackster_eta[valType].push_back(ibook.book1D(
1359 "NumDup_Trackster_Eta" + valSuffix_[valType], "Num Duplicate Trackster vs Eta;#eta", nintEta_, minEta_, maxEta_));
1360 histograms.h_denom_caloparticle_eta[valType].push_back(
1361 ibook.book1D("Denom_" + ref_[valType] + "_Eta",
1362 "Denom " + refText_[valType] + " Eta per Trackster;#eta",
1363 nintEta_,
1364 minEta_,
1365 maxEta_));
1366
1367 histograms.h_numEff_caloparticle_phi[valType].push_back(
1368 ibook.book1D("NumEff_" + ref_[valType] + "_Phi",
1369 "Num Efficiency " + refText_[valType] + " Phi per Trackster;#phi",
1370 nintPhi_,
1371 minPhi_,
1372 maxPhi_));
1373 histograms.h_num_caloparticle_phi[valType].push_back(
1374 ibook.book1D("Num_" + ref_[valType] + "_Phi",
1375 "Num Purity " + refText_[valType] + " Phi per Trackster;#phi",
1376 nintPhi_,
1377 minPhi_,
1378 maxPhi_));
1379 histograms.h_numDup_trackster_phi[valType].push_back(ibook.book1D(
1380 "NumDup_Trackster_Phi" + valSuffix_[valType], "Num Duplicate Trackster vs Phi;#phi", nintPhi_, minPhi_, maxPhi_));
1381 histograms.h_denom_caloparticle_phi[valType].push_back(
1382 ibook.book1D("Denom_" + ref_[valType] + "_Phi",
1383 "Denom " + refText_[valType] + " Phi per Trackster;#phi",
1384 nintPhi_,
1385 minPhi_,
1386 maxPhi_));
1387
1388 histograms.h_numEff_caloparticle_en[valType].push_back(
1389 ibook.book1D("NumEff_" + ref_[valType] + "_Energy",
1390 "Num Efficiency " + refText_[valType] + " Energy per Trackster;energy [GeV]",
1391 nintEne_,
1392 minEne_,
1393 maxEne_));
1394 histograms.h_num_caloparticle_en[valType].push_back(
1395 ibook.book1D("Num_" + ref_[valType] + "_Energy",
1396 "Num Purity " + refText_[valType] + " Energy per Trackster;energy [GeV]",
1397 nintEne_,
1398 minEne_,
1399 maxEne_));
1400 histograms.h_numDup_trackster_en[valType].push_back(ibook.book1D("NumDup_Trackster_Energy" + valSuffix_[valType],
1401 "Num Duplicate Trackster vs Energy;energy [GeV]",
1402 nintEne_,
1403 minEne_,
1404 maxEne_));
1405 histograms.h_denom_caloparticle_en[valType].push_back(
1406 ibook.book1D("Denom_" + ref_[valType] + "_Energy",
1407 "Denom " + refText_[valType] + " Energy per Trackster;energy [GeV]",
1408 nintEne_,
1409 minEne_,
1410 maxEne_));
1411
1412 histograms.h_numEff_caloparticle_pt[valType].push_back(
1413 ibook.book1D("NumEff_" + ref_[valType] + "_Pt",
1414 "Num Efficiency " + refText_[valType] + " p_{T} per Trackster;p_{T} [GeV]",
1415 nintPt_,
1416 minPt_,
1417 maxPt_));
1418 histograms.h_num_caloparticle_pt[valType].push_back(
1419 ibook.book1D("Num_" + ref_[valType] + "_Pt",
1420 "Num Purity " + refText_[valType] + " p_{T} per Trackster;p_{T} [GeV]",
1421 nintPt_,
1422 minPt_,
1423 maxPt_));
1424 histograms.h_numDup_trackster_pt[valType].push_back(ibook.book1D("NumDup_Trackster_Pt" + valSuffix_[valType],
1425 "Num Duplicate Trackster vs p_{T};p_{T} [GeV]",
1426 nintPt_,
1427 minPt_,
1428 maxPt_));
1429 histograms.h_denom_caloparticle_pt[valType].push_back(
1430 ibook.book1D("Denom_" + ref_[valType] + "_Pt",
1431 "Denom " + refText_[valType] + " p_{T} per Trackster;p_{T} [GeV]",
1432 nintPt_,
1433 minPt_,
1434 maxPt_));
1435 }
1436
1437 void HGVHistoProducerAlgo::fill_info_histos(const Histograms& histograms, unsigned int layers) const {
1438
1439
1440
1441
1442 histograms.lastLayerEEzm->Fill(recHitTools_->lastLayerEE());
1443 histograms.lastLayerFHzm->Fill(recHitTools_->lastLayerFH());
1444 histograms.maxlayerzm->Fill(layers);
1445 histograms.lastLayerEEzp->Fill(recHitTools_->lastLayerEE() + layers);
1446 histograms.lastLayerFHzp->Fill(recHitTools_->lastLayerFH() + layers);
1447 histograms.maxlayerzp->Fill(layers + layers);
1448 }
1449
1450 void HGVHistoProducerAlgo::fill_caloparticle_histos(const Histograms& histograms,
1451 int pdgid,
1452 const CaloParticle& caloParticle,
1453 std::vector<SimVertex> const& simVertices,
1454 unsigned int layers,
1455 std::unordered_map<DetId, const unsigned int> const& hitMap,
1456 MultiVectorManager<HGCRecHit> const& hits) const {
1457 const auto eta = getEta(caloParticle.eta());
1458 if (histograms.h_caloparticle_eta.count(pdgid)) {
1459 histograms.h_caloparticle_eta.at(pdgid)->Fill(eta);
1460 }
1461 if (histograms.h_caloparticle_eta_Zorigin.count(pdgid)) {
1462 histograms.h_caloparticle_eta_Zorigin.at(pdgid)->Fill(
1463 simVertices.at(caloParticle.g4Tracks()[0].vertIndex()).position().z(), eta);
1464 }
1465
1466 if (histograms.h_caloparticle_energy.count(pdgid)) {
1467 histograms.h_caloparticle_energy.at(pdgid)->Fill(caloParticle.energy());
1468 }
1469 if (histograms.h_caloparticle_pt.count(pdgid)) {
1470 histograms.h_caloparticle_pt.at(pdgid)->Fill(caloParticle.pt());
1471 }
1472 if (histograms.h_caloparticle_phi.count(pdgid)) {
1473 histograms.h_caloparticle_phi.at(pdgid)->Fill(caloParticle.phi());
1474 }
1475
1476 if (histograms.h_caloparticle_nSimClusters.count(pdgid)) {
1477 histograms.h_caloparticle_nSimClusters.at(pdgid)->Fill(caloParticle.simClusters().size());
1478
1479 int simHits = 0;
1480 int minLayerId = 999;
1481 int maxLayerId = 0;
1482
1483 int simHits_matched = 0;
1484 int minLayerId_matched = 999;
1485 int maxLayerId_matched = 0;
1486
1487 float energy = 0.;
1488 std::map<int, double> totenergy_layer;
1489
1490 float hitEnergyWeight_invSum = 0;
1491 std::vector<std::pair<DetId, float>> haf_cp;
1492 for (const auto& sc : caloParticle.simClusters()) {
1493 LogDebug("HGCalValidator") << " This sim cluster has " << sc->hits_and_fractions().size() << " simHits and "
1494 << sc->energy() << " energy. " << std::endl;
1495 simHits += sc->endcap_hits_and_fractions().size();
1496 for (auto const& h_and_f : sc->endcap_hits_and_fractions()) {
1497 const auto hitDetId = h_and_f.first;
1498 if (recHitTools_->isBarrel(hitDetId))
1499 continue;
1500 const int layerId =
1501 recHitTools_->getLayerWithOffset(hitDetId) + layers * ((recHitTools_->zside(hitDetId) + 1) >> 1) - 1;
1502
1503 int layerId_matched_min = 999;
1504 int layerId_matched_max = 0;
1505 std::unordered_map<DetId, const unsigned int>::const_iterator itcheck = hitMap.find(hitDetId);
1506 if (itcheck != hitMap.end()) {
1507 layerId_matched_min = layerId;
1508 layerId_matched_max = layerId;
1509 simHits_matched++;
1510
1511 const auto hitEn = (hits[itcheck->second]).energy();
1512 hitEnergyWeight_invSum += pow(hitEn, 2);
1513 const auto hitFr = h_and_f.second;
1514 const auto hitEnFr = hitEn * hitFr;
1515 energy += hitEnFr;
1516 histograms.h_caloparticle_nHits_matched_energy.at(pdgid)->Fill(hitEnFr);
1517 histograms.h_caloparticle_nHits_matched_energy_layer.at(pdgid)->Fill(layerId, hitEnFr);
1518
1519 if (totenergy_layer.find(layerId) != totenergy_layer.end()) {
1520 totenergy_layer[layerId] = totenergy_layer.at(layerId) + hitEn;
1521 } else {
1522 totenergy_layer.emplace(layerId, hitEn);
1523 }
1524 if (caloParticle.simClusters().size() == 1)
1525 histograms.h_caloparticle_nHits_matched_energy_layer_1SimCl.at(pdgid)->Fill(layerId, hitEnFr);
1526
1527 auto found = std::find_if(std::begin(haf_cp),
1528 std::end(haf_cp),
1529 [&hitDetId](const std::pair<DetId, float>& v) { return v.first == hitDetId; });
1530 if (found != haf_cp.end())
1531 found->second += hitFr;
1532 else
1533 haf_cp.emplace_back(hitDetId, hitFr);
1534
1535 } else {
1536 LogDebug("HGCalValidator") << " matched to RecHit NOT found !" << std::endl;
1537 }
1538
1539 minLayerId = std::min(minLayerId, layerId);
1540 maxLayerId = std::max(maxLayerId, layerId);
1541 minLayerId_matched = std::min(minLayerId_matched, layerId_matched_min);
1542 maxLayerId_matched = std::max(maxLayerId_matched, layerId_matched_max);
1543 }
1544 LogDebug("HGCalValidator") << std::endl;
1545 }
1546 if (hitEnergyWeight_invSum)
1547 hitEnergyWeight_invSum = 1 / hitEnergyWeight_invSum;
1548
1549 if (minLayerId == 999)
1550 return;
1551 histograms.h_caloparticle_firstlayer.at(pdgid)->Fill(minLayerId);
1552 histograms.h_caloparticle_lastlayer.at(pdgid)->Fill(maxLayerId);
1553 histograms.h_caloparticle_layersnum.at(pdgid)->Fill(int(maxLayerId - minLayerId));
1554
1555 histograms.h_caloparticle_firstlayer_matchedtoRecHit.at(pdgid)->Fill(minLayerId_matched);
1556 histograms.h_caloparticle_lastlayer_matchedtoRecHit.at(pdgid)->Fill(maxLayerId_matched);
1557 histograms.h_caloparticle_layersnum_matchedtoRecHit.at(pdgid)->Fill(int(maxLayerId_matched - minLayerId_matched));
1558
1559 histograms.h_caloparticle_nHitsInSimClusters.at(pdgid)->Fill((float)simHits);
1560 histograms.h_caloparticle_nHitsInSimClusters_matchedtoRecHit.at(pdgid)->Fill((float)simHits_matched);
1561 histograms.h_caloparticle_selfenergy.at(pdgid)->Fill((float)energy);
1562 histograms.h_caloparticle_energyDifference.at(pdgid)->Fill((float)1. - energy / caloParticle.energy());
1563
1564
1565 auto i = totenergy_layer.begin();
1566 double sum_energy = 0.0;
1567 while (i != totenergy_layer.end()) {
1568 sum_energy += i->second;
1569 histograms.h_caloparticle_sum_energy_layer.at(pdgid)->Fill(i->first, sum_energy / caloParticle.energy() * 100.);
1570 i++;
1571 }
1572
1573 for (auto const& haf : haf_cp) {
1574 const auto hitEn = (hits[hitMap.find(haf.first)->second]).energy();
1575 const auto weight = pow(hitEn, 2);
1576 histograms.h_caloparticle_fractions.at(pdgid)->Fill(haf.second, weight * hitEnergyWeight_invSum);
1577 histograms.h_caloparticle_fractions_weight.at(pdgid)->Fill(haf.second, weight * hitEnergyWeight_invSum, weight);
1578 }
1579 }
1580 }
1581
1582 void HGVHistoProducerAlgo::HGVHistoProducerAlgo::fill_simCluster_histos(const Histograms& histograms,
1583 std::vector<SimCluster> const& simClusters,
1584 unsigned int layers,
1585 std::vector<int> thicknesses) const {
1586
1587
1588
1589
1590
1591
1592
1593
1594 std::vector<int> tnscpl(1000, 0);
1595
1596
1597 std::map<std::string, int> tnscpthplus;
1598 tnscpthplus.clear();
1599 std::map<std::string, int> tnscpthminus;
1600 tnscpthminus.clear();
1601
1602 for (std::vector<int>::iterator it = thicknesses.begin(); it != thicknesses.end(); ++it) {
1603 tnscpthplus.insert(std::pair<std::string, int>(std::to_string(*it), 0));
1604 tnscpthminus.insert(std::pair<std::string, int>(std::to_string(*it), 0));
1605 }
1606
1607 tnscpthplus.insert(std::pair<std::string, int>("mixed", 0));
1608 tnscpthminus.insert(std::pair<std::string, int>("mixed", 0));
1609
1610
1611 for (const auto& sc : simClusters) {
1612
1613 int nthhits120p = 0;
1614 int nthhits200p = 0;
1615 int nthhits300p = 0;
1616 int nthhitsscintp = 0;
1617 int nthhits120m = 0;
1618 int nthhits200m = 0;
1619 int nthhits300m = 0;
1620 int nthhitsscintm = 0;
1621
1622 double thickness = 0.;
1623
1624 std::vector<int> occurenceSCinlayer(1000, 0);
1625
1626
1627 for (const auto& hAndF : sc.hits_and_fractions()) {
1628 const DetId sh_detid = hAndF.first;
1629
1630 if (sh_detid.det() == DetId::Forward || sh_detid.det() == DetId::HGCalEE || sh_detid.det() == DetId::HGCalHSi ||
1631 sh_detid.det() == DetId::HGCalHSc) {
1632
1633 int layerid =
1634 recHitTools_->getLayerWithOffset(sh_detid) + layers * ((recHitTools_->zside(sh_detid) + 1) >> 1) - 1;
1635
1636 int zside = recHitTools_->zside(sh_detid);
1637
1638
1639 if (occurenceSCinlayer[layerid] == 0) {
1640 tnscpl[layerid]++;
1641 }
1642 occurenceSCinlayer[layerid]++;
1643
1644 if (sh_detid.det() == DetId::HGCalHSc)
1645 thickness = -1;
1646 else
1647 thickness = recHitTools_->getSiThickness(sh_detid);
1648
1649 if ((thickness == 120.) && (zside > 0.)) {
1650 nthhits120p++;
1651 } else if ((thickness == 120.) && (zside < 0.)) {
1652 nthhits120m++;
1653 } else if ((thickness == 200.) && (zside > 0.)) {
1654 nthhits200p++;
1655 } else if ((thickness == 200.) && (zside < 0.)) {
1656 nthhits200m++;
1657 } else if ((thickness == 300.) && (zside > 0.)) {
1658 nthhits300p++;
1659 } else if ((thickness == 300.) && (zside < 0.)) {
1660 nthhits300m++;
1661 } else if ((thickness == -1) && (zside > 0.)) {
1662 nthhitsscintp++;
1663 } else if ((thickness == -1) && (zside < 0.)) {
1664 nthhitsscintm++;
1665 } else {
1666 LogDebug("HGCalValidator")
1667 << " You are running a geometry that contains thicknesses different than the normal ones. "
1668 << "\n";
1669 }
1670 }
1671 }
1672
1673
1674 if ((nthhits120p != 0 && nthhits200p != 0) || (nthhits120p != 0 && nthhits300p != 0) ||
1675 (nthhits120p != 0 && nthhitsscintp != 0) || (nthhits200p != 0 && nthhits300p != 0) ||
1676 (nthhits200p != 0 && nthhitsscintp != 0) || (nthhits300p != 0 && nthhitsscintp != 0)) {
1677 tnscpthplus["mixed"]++;
1678 } else if ((nthhits120p != 0 || nthhits200p != 0 || nthhits300p != 0 || nthhitsscintp != 0)) {
1679
1680 tnscpthplus[std::to_string((int)thickness)]++;
1681 }
1682 if ((nthhits120m != 0 && nthhits200m != 0) || (nthhits120m != 0 && nthhits300m != 0) ||
1683 (nthhits120m != 0 && nthhitsscintm != 0) || (nthhits200m != 0 && nthhits300m != 0) ||
1684 (nthhits200m != 0 && nthhitsscintm != 0) || (nthhits300m != 0 && nthhitsscintm != 0)) {
1685 tnscpthminus["mixed"]++;
1686 } else if ((nthhits120m != 0 || nthhits200m != 0 || nthhits300m != 0 || nthhitsscintm != 0)) {
1687
1688 tnscpthminus[std::to_string((int)thickness)]++;
1689 }
1690 }
1691
1692
1693 for (unsigned ilayer = 0; ilayer < layers * 2; ++ilayer)
1694 if (histograms.h_simclusternum_perlayer.count(ilayer))
1695 histograms.h_simclusternum_perlayer.at(ilayer)->Fill(tnscpl[ilayer]);
1696
1697
1698 for (std::vector<int>::iterator it = thicknesses.begin(); it != thicknesses.end(); ++it) {
1699 if (histograms.h_simclusternum_perthick.count(*it)) {
1700 histograms.h_simclusternum_perthick.at(*it)->Fill(tnscpthplus[std::to_string(*it)]);
1701 histograms.h_simclusternum_perthick.at(*it)->Fill(tnscpthminus[std::to_string(*it)]);
1702 }
1703 }
1704
1705 histograms.h_mixedhitssimcluster_zplus->Fill(tnscpthplus["mixed"]);
1706 histograms.h_mixedhitssimcluster_zminus->Fill(tnscpthminus["mixed"]);
1707 }
1708
1709 void HGVHistoProducerAlgo::HGVHistoProducerAlgo::fill_simClusterAssociation_histos(
1710 const Histograms& histograms,
1711 const int count,
1712 edm::Handle<reco::CaloClusterCollection> clusterHandle,
1713 const reco::CaloClusterCollection& clusters,
1714 edm::Handle<std::vector<SimCluster>> simClusterHandle,
1715 std::vector<SimCluster> const& simClusters,
1716 std::vector<size_t> const& sCIndices,
1717 const std::vector<float>& mask,
1718 std::unordered_map<DetId, const unsigned int> const& hitMap,
1719 unsigned int layers,
1720 const ticl::RecoToSimCollectionWithSimClusters& scsInLayerClusterMap,
1721 const ticl::SimToRecoCollectionWithSimClusters& lcsInSimClusterMap,
1722 MultiVectorManager<HGCRecHit> const& hits) const {
1723
1724
1725
1726
1727
1728
1729
1730
1731 layerClusters_to_SimClusters(histograms,
1732 count,
1733 clusterHandle,
1734 clusters,
1735 simClusterHandle,
1736 simClusters,
1737 sCIndices,
1738 mask,
1739 hitMap,
1740 layers,
1741 scsInLayerClusterMap,
1742 lcsInSimClusterMap,
1743 hits);
1744 }
1745
1746 void HGVHistoProducerAlgo::fill_cluster_histos(const Histograms& histograms,
1747 const int count,
1748 const reco::CaloCluster& cluster) const {
1749 const auto eta = getEta(cluster.eta());
1750 histograms.h_cluster_eta[count]->Fill(eta);
1751 }
1752
1753 void HGVHistoProducerAlgo::layerClusters_to_CaloParticles(const Histograms& histograms,
1754 edm::Handle<reco::CaloClusterCollection> clusterHandle,
1755 const reco::CaloClusterCollection& clusters,
1756 edm::Handle<std::vector<CaloParticle>> caloParticleHandle,
1757 std::vector<CaloParticle> const& cP,
1758 std::vector<size_t> const& cPIndices,
1759 std::vector<size_t> const& cPSelectedIndices,
1760 std::unordered_map<DetId, const unsigned int> const& hitMap,
1761 unsigned int layers,
1762 const ticl::RecoToSimCollection& cpsInLayerClusterMap,
1763 const ticl::SimToRecoCollection& cPOnLayerMap,
1764 MultiVectorManager<HGCRecHit> const& hits) const {
1765 const auto nLayerClusters = clusters.size();
1766
1767 std::unordered_map<DetId, std::vector<HGVHistoProducerAlgo::detIdInfoInCluster>> detIdToCaloParticleId_Map;
1768 std::unordered_map<DetId, std::vector<HGVHistoProducerAlgo::detIdInfoInCluster>> detIdToLayerClusterId_Map;
1769
1770
1771
1772 for (const auto& cpId : cPIndices) {
1773 for (const auto& simCluster : cP[cpId].simClusters()) {
1774 for (const auto& it_haf : simCluster->hits_and_fractions()) {
1775 const DetId hitid = (it_haf.first);
1776 if (recHitTools_->isBarrel(hitid))
1777 continue;
1778 if (hitMap.find(hitid) != hitMap.end()) {
1779 if (detIdToCaloParticleId_Map.find(hitid) == detIdToCaloParticleId_Map.end()) {
1780 detIdToCaloParticleId_Map[hitid] = std::vector<HGVHistoProducerAlgo::detIdInfoInCluster>();
1781 detIdToCaloParticleId_Map[hitid].emplace_back(
1782 HGVHistoProducerAlgo::detIdInfoInCluster{cpId, it_haf.second});
1783 } else {
1784 auto findHitIt =
1785 std::find(detIdToCaloParticleId_Map[hitid].begin(),
1786 detIdToCaloParticleId_Map[hitid].end(),
1787 HGVHistoProducerAlgo::detIdInfoInCluster{
1788 cpId, 0.f});
1789 if (findHitIt != detIdToCaloParticleId_Map[hitid].end())
1790 findHitIt->fraction += it_haf.second;
1791 else
1792 detIdToCaloParticleId_Map[hitid].emplace_back(
1793 HGVHistoProducerAlgo::detIdInfoInCluster{cpId, it_haf.second});
1794 }
1795 }
1796 }
1797 }
1798 }
1799
1800 for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) {
1801 const auto& hits_and_fractions = clusters[lcId].hitsAndFractions();
1802 const auto numberOfHitsInLC = hits_and_fractions.size();
1803
1804
1805
1806
1807
1808
1809
1810
1811
1812 std::vector<int> hitsToCaloParticleId(numberOfHitsInLC);
1813 const auto firstHitDetId = hits_and_fractions[0].first;
1814 if (recHitTools_->isBarrel(firstHitDetId))
1815 continue;
1816 int lcLayerId =
1817 recHitTools_->getLayerWithOffset(firstHitDetId) + layers * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
1818
1819
1820 std::unordered_map<unsigned, float> CPEnergyInLC;
1821
1822 for (unsigned int iHit = 0; iHit < numberOfHitsInLC; iHit++) {
1823 const DetId rh_detid = hits_and_fractions[iHit].first;
1824 const auto rhFraction = hits_and_fractions[iHit].second;
1825
1826 std::unordered_map<DetId, const unsigned int>::const_iterator itcheck = hitMap.find(rh_detid);
1827 const HGCRecHit* hit = &(hits[itcheck->second]);
1828
1829 if (detIdToLayerClusterId_Map.find(rh_detid) == detIdToLayerClusterId_Map.end()) {
1830 detIdToLayerClusterId_Map[rh_detid] = std::vector<HGVHistoProducerAlgo::detIdInfoInCluster>();
1831 }
1832 detIdToLayerClusterId_Map[rh_detid].emplace_back(HGVHistoProducerAlgo::detIdInfoInCluster{lcId, rhFraction});
1833
1834 const auto& hit_find_in_CP = detIdToCaloParticleId_Map.find(rh_detid);
1835
1836
1837
1838
1839
1840
1841
1842 if (rhFraction == 0.) {
1843 hitsToCaloParticleId[iHit] = -2;
1844 }
1845 if (hit_find_in_CP == detIdToCaloParticleId_Map.end()) {
1846 hitsToCaloParticleId[iHit] -= 1;
1847 } else {
1848 auto maxCPEnergyInLC = 0.f;
1849 auto maxCPId = -1;
1850 for (auto& h : hit_find_in_CP->second) {
1851 const auto iCP = h.clusterId;
1852 CPEnergyInLC[iCP] += h.fraction * hit->energy();
1853
1854
1855 if (CPEnergyInLC[iCP] > maxCPEnergyInLC) {
1856 maxCPEnergyInLC = CPEnergyInLC[iCP];
1857 maxCPId = iCP;
1858 }
1859 }
1860 hitsToCaloParticleId[iHit] = maxCPId;
1861 }
1862 histograms.h_cellAssociation_perlayer.at(lcLayerId)->Fill(
1863 hitsToCaloParticleId[iHit] > 0. ? 0. : hitsToCaloParticleId[iHit]);
1864 }
1865
1866 }
1867
1868
1869
1870
1871 for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) {
1872 const auto firstHitDetId = (clusters[lcId].hitsAndFractions())[0].first;
1873 if (recHitTools_->isBarrel(firstHitDetId))
1874 continue;
1875 const int lcLayerId =
1876 recHitTools_->getLayerWithOffset(firstHitDetId) + layers * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
1877 histograms.h_denom_layercl_eta_perlayer.at(lcLayerId)->Fill(clusters[lcId].eta());
1878 histograms.h_denom_layercl_phi_perlayer.at(lcLayerId)->Fill(clusters[lcId].phi());
1879
1880 const edm::Ref<reco::CaloClusterCollection> lcRef(clusterHandle, lcId);
1881 const auto& cpsIt = cpsInLayerClusterMap.find(lcRef);
1882 if (cpsIt == cpsInLayerClusterMap.end())
1883 continue;
1884
1885 const auto lc_en = clusters[lcId].energy();
1886 const auto& cps = cpsIt->val;
1887 if (lc_en == 0. && !cps.empty()) {
1888 for (const auto& cpPair : cps)
1889 histograms.h_score_layercl2caloparticle_perlayer.at(lcLayerId)->Fill(cpPair.second);
1890 continue;
1891 }
1892 for (const auto& cpPair : cps) {
1893 LogDebug("HGCalValidator") << "layerCluster Id: \t" << lcId << "\t CP id: \t" << cpPair.first.index()
1894 << "\t score \t" << cpPair.second << std::endl;
1895 histograms.h_score_layercl2caloparticle_perlayer.at(lcLayerId)->Fill(cpPair.second);
1896 auto const& cp_linked =
1897 std::find_if(std::begin(cPOnLayerMap[cpPair.first]),
1898 std::end(cPOnLayerMap[cpPair.first]),
1899 [&lcRef](const std::pair<edm::Ref<reco::CaloClusterCollection>, std::pair<float, float>>& p) {
1900 return p.first == lcRef;
1901 });
1902 if (cp_linked ==
1903 cPOnLayerMap[cpPair.first].end())
1904 continue;
1905 histograms.h_sharedenergy_layercl2caloparticle_perlayer.at(lcLayerId)->Fill(cp_linked->second.first / lc_en,
1906 lc_en);
1907 histograms.h_energy_vs_score_layercl2caloparticle_perlayer.at(lcLayerId)->Fill(cpPair.second,
1908 cp_linked->second.first / lc_en);
1909 }
1910 const auto assoc =
1911 std::count_if(std::begin(cps), std::end(cps), [](const auto& obj) { return obj.second < ScoreCutLCtoCP_; });
1912 if (assoc) {
1913 histograms.h_num_layercl_eta_perlayer.at(lcLayerId)->Fill(clusters[lcId].eta());
1914 histograms.h_num_layercl_phi_perlayer.at(lcLayerId)->Fill(clusters[lcId].phi());
1915 if (assoc > 1) {
1916 histograms.h_numMerge_layercl_eta_perlayer.at(lcLayerId)->Fill(clusters[lcId].eta());
1917 histograms.h_numMerge_layercl_phi_perlayer.at(lcLayerId)->Fill(clusters[lcId].phi());
1918 }
1919 const auto& best = std::min_element(
1920 std::begin(cps), std::end(cps), [](const auto& obj1, const auto& obj2) { return obj1.second < obj2.second; });
1921 const auto& best_cp_linked =
1922 std::find_if(std::begin(cPOnLayerMap[best->first]),
1923 std::end(cPOnLayerMap[best->first]),
1924 [&lcRef](const std::pair<edm::Ref<reco::CaloClusterCollection>, std::pair<float, float>>& p) {
1925 return p.first == lcRef;
1926 });
1927 if (best_cp_linked ==
1928 cPOnLayerMap[best->first].end())
1929 continue;
1930 histograms.h_sharedenergy_layercl2caloparticle_vs_eta_perlayer.at(lcLayerId)->Fill(
1931 clusters[lcId].eta(), best_cp_linked->second.first / lc_en);
1932 histograms.h_sharedenergy_layercl2caloparticle_vs_phi_perlayer.at(lcLayerId)->Fill(
1933 clusters[lcId].phi(), best_cp_linked->second.first / lc_en);
1934 }
1935 }
1936
1937
1938
1939
1940 for (const auto& cpId : cPSelectedIndices) {
1941 const edm::Ref<CaloParticleCollection> cpRef(caloParticleHandle, cpId);
1942 const auto& lcsIt = cPOnLayerMap.find(cpRef);
1943
1944 std::map<unsigned int, float> cPEnergyOnLayer;
1945 for (unsigned int layerId = 0; layerId < layers * 2; ++layerId)
1946 cPEnergyOnLayer[layerId] = 0;
1947
1948 for (const auto& simCluster : cP[cpId].simClusters()) {
1949 for (const auto& it_haf : simCluster->hits_and_fractions()) {
1950 const DetId hitid = (it_haf.first);
1951 if (recHitTools_->isBarrel(hitid))
1952 continue;
1953 const auto hitLayerId =
1954 recHitTools_->getLayerWithOffset(hitid) + layers * ((recHitTools_->zside(hitid) + 1) >> 1) - 1;
1955 std::unordered_map<DetId, const unsigned int>::const_iterator itcheck = hitMap.find(hitid);
1956 if (itcheck != hitMap.end()) {
1957 const HGCRecHit* hit = &(hits[itcheck->second]);
1958 cPEnergyOnLayer[hitLayerId] += it_haf.second * hit->energy();
1959 }
1960 }
1961 }
1962
1963 for (unsigned int layerId = 0; layerId < layers * 2; ++layerId) {
1964 if (!cPEnergyOnLayer[layerId])
1965 continue;
1966
1967 histograms.h_denom_caloparticle_eta_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().eta());
1968 histograms.h_denom_caloparticle_phi_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().phi());
1969
1970 if (lcsIt == cPOnLayerMap.end())
1971 continue;
1972 const auto& lcs = lcsIt->val;
1973
1974 auto getLCLayerId = [&](const unsigned int lcId) {
1975 const auto firstHitDetId = (clusters[lcId].hitsAndFractions())[0].first;
1976 if (recHitTools_->isBarrel(firstHitDetId))
1977 return 9999u;
1978 const auto lcLayerId = recHitTools_->getLayerWithOffset(firstHitDetId) +
1979 layers * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
1980 return lcLayerId;
1981 };
1982
1983 for (const auto& lcPair : lcs) {
1984 if (recHitTools_->isBarrel(clusters[lcPair.first.index()].seed()))
1985 continue;
1986 if (getLCLayerId(lcPair.first.index()) != layerId)
1987 continue;
1988 histograms.h_score_caloparticle2layercl_perlayer.at(layerId)->Fill(lcPair.second.second);
1989 histograms.h_sharedenergy_caloparticle2layercl_perlayer.at(layerId)->Fill(
1990 lcPair.second.first / cPEnergyOnLayer[layerId], cPEnergyOnLayer[layerId]);
1991 histograms.h_energy_vs_score_caloparticle2layercl_perlayer.at(layerId)->Fill(
1992 lcPair.second.second, lcPair.second.first / cPEnergyOnLayer[layerId]);
1993 }
1994 const auto assoc = std::count_if(std::begin(lcs), std::end(lcs), [&](const auto& obj) {
1995 if (getLCLayerId(obj.first.index()) != layerId)
1996 return false;
1997 else
1998 return obj.second.second < ScoreCutCPtoLC_;
1999 });
2000 if (assoc) {
2001 histograms.h_num_caloparticle_eta_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().eta());
2002 histograms.h_num_caloparticle_phi_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().phi());
2003 if (assoc > 1) {
2004 histograms.h_numDup_caloparticle_eta_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().eta());
2005 histograms.h_numDup_caloparticle_phi_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().phi());
2006 }
2007 const auto best = std::min_element(std::begin(lcs), std::end(lcs), [&](const auto& obj1, const auto& obj2) {
2008 if (getLCLayerId(obj1.first.index()) != layerId)
2009 return false;
2010 else if (getLCLayerId(obj2.first.index()) == layerId)
2011 return obj1.second.second < obj2.second.second;
2012 else
2013 return true;
2014 });
2015 histograms.h_sharedenergy_caloparticle2layercl_vs_eta_perlayer.at(layerId)->Fill(
2016 cP[cpId].g4Tracks()[0].momentum().eta(), best->second.first / cPEnergyOnLayer[layerId]);
2017 histograms.h_sharedenergy_caloparticle2layercl_vs_phi_perlayer.at(layerId)->Fill(
2018 cP[cpId].g4Tracks()[0].momentum().phi(), best->second.first / cPEnergyOnLayer[layerId]);
2019 }
2020 }
2021 }
2022 }
2023
2024 void HGVHistoProducerAlgo::layerClusters_to_SimClusters(
2025 const Histograms& histograms,
2026 const int count,
2027 edm::Handle<reco::CaloClusterCollection> clusterHandle,
2028 const reco::CaloClusterCollection& clusters,
2029 edm::Handle<std::vector<SimCluster>> simClusterHandle,
2030 std::vector<SimCluster> const& sC,
2031 std::vector<size_t> const& sCIndices,
2032 const std::vector<float>& mask,
2033 std::unordered_map<DetId, const unsigned int> const& hitMap,
2034 unsigned int layers,
2035 const ticl::RecoToSimCollectionWithSimClusters& scsInLayerClusterMap,
2036 const ticl::SimToRecoCollectionWithSimClusters& lcsInSimClusterMap,
2037 MultiVectorManager<HGCRecHit> const& hits) const {
2038
2039
2040
2041 for (unsigned int lcId = 0; lcId < clusters.size(); ++lcId) {
2042 if (mask[lcId] != 0.) {
2043 LogDebug("HGCalValidator") << "Skipping layer cluster " << lcId << " not belonging to mask" << std::endl;
2044 continue;
2045 }
2046 const auto firstHitDetId = (clusters[lcId].hitsAndFractions())[0].first;
2047 if (recHitTools_->isBarrel(firstHitDetId))
2048 continue;
2049 const auto lcLayerId =
2050 recHitTools_->getLayerWithOffset(firstHitDetId) + layers * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
2051
2052
2053 histograms.h_denom_layercl_in_simcl_eta_perlayer[count].at(lcLayerId)->Fill(clusters[lcId].eta());
2054 histograms.h_denom_layercl_in_simcl_phi_perlayer[count].at(lcLayerId)->Fill(clusters[lcId].phi());
2055
2056 const edm::Ref<reco::CaloClusterCollection> lcRef(clusterHandle, lcId);
2057 const auto& scsIt = scsInLayerClusterMap.find(lcRef);
2058 if (scsIt == scsInLayerClusterMap.end())
2059 continue;
2060
2061 const auto lc_en = clusters[lcId].energy();
2062 const auto& scs = scsIt->val;
2063
2064
2065 if (lc_en == 0. && !scs.empty()) {
2066 for (const auto& scPair : scs) {
2067 histograms.h_score_layercl2simcluster_perlayer[count].at(lcLayerId)->Fill(scPair.second);
2068 }
2069 continue;
2070 }
2071
2072 for (const auto& scPair : scs) {
2073 LogDebug("HGCalValidator") << "layerCluster Id: \t" << lcId << "\t SC id: \t" << scPair.first.index()
2074 << "\t score \t" << scPair.second << std::endl;
2075
2076 histograms.h_score_layercl2simcluster_perlayer[count].at(lcLayerId)->Fill(scPair.second);
2077 auto const& sc_linked =
2078 std::find_if(std::begin(lcsInSimClusterMap[scPair.first]),
2079 std::end(lcsInSimClusterMap[scPair.first]),
2080 [&lcRef](const std::pair<edm::Ref<reco::CaloClusterCollection>, std::pair<float, float>>& p) {
2081 return p.first == lcRef;
2082 });
2083 if (sc_linked ==
2084 lcsInSimClusterMap[scPair.first].end())
2085 continue;
2086 histograms.h_sharedenergy_layercl2simcluster_perlayer[count].at(lcLayerId)->Fill(sc_linked->second.first / lc_en,
2087 lc_en);
2088 histograms.h_energy_vs_score_layercl2simcluster_perlayer[count].at(lcLayerId)->Fill(
2089 scPair.second, sc_linked->second.first / lc_en);
2090 }
2091
2092 const auto assoc =
2093 std::count_if(std::begin(scs), std::end(scs), [](const auto& obj) { return obj.second < ScoreCutLCtoSC_; });
2094 if (assoc) {
2095 histograms.h_num_layercl_in_simcl_eta_perlayer[count].at(lcLayerId)->Fill(clusters[lcId].eta());
2096 histograms.h_num_layercl_in_simcl_phi_perlayer[count].at(lcLayerId)->Fill(clusters[lcId].phi());
2097 if (assoc > 1) {
2098 histograms.h_numMerge_layercl_in_simcl_eta_perlayer[count].at(lcLayerId)->Fill(clusters[lcId].eta());
2099 histograms.h_numMerge_layercl_in_simcl_phi_perlayer[count].at(lcLayerId)->Fill(clusters[lcId].phi());
2100 }
2101 const auto& best = std::min_element(
2102 std::begin(scs), std::end(scs), [](const auto& obj1, const auto& obj2) { return obj1.second < obj2.second; });
2103
2104 const auto& best_sc_linked =
2105 std::find_if(std::begin(lcsInSimClusterMap[best->first]),
2106 std::end(lcsInSimClusterMap[best->first]),
2107 [&lcRef](const std::pair<edm::Ref<reco::CaloClusterCollection>, std::pair<float, float>>& p) {
2108 return p.first == lcRef;
2109 });
2110 if (best_sc_linked ==
2111 lcsInSimClusterMap[best->first].end())
2112 continue;
2113 histograms.h_sharedenergy_layercl2simcluster_vs_eta_perlayer[count].at(lcLayerId)->Fill(
2114 clusters[lcId].eta(), best_sc_linked->second.first / lc_en);
2115 histograms.h_sharedenergy_layercl2simcluster_vs_phi_perlayer[count].at(lcLayerId)->Fill(
2116 clusters[lcId].phi(), best_sc_linked->second.first / lc_en);
2117 }
2118 }
2119
2120
2121
2122
2123 for (const auto& scId : sCIndices) {
2124 const edm::Ref<SimClusterCollection> scRef(simClusterHandle, scId);
2125 const auto& lcsIt = lcsInSimClusterMap.find(scRef);
2126
2127 std::map<unsigned int, float> sCEnergyOnLayer;
2128 for (unsigned int layerId = 0; layerId < layers * 2; ++layerId)
2129 sCEnergyOnLayer[layerId] = 0;
2130
2131 for (const auto& it_haf : sC[scId].hits_and_fractions()) {
2132 const DetId hitid = (it_haf.first);
2133 if (recHitTools_->isBarrel(hitid))
2134 continue;
2135 const auto scLayerId =
2136 recHitTools_->getLayerWithOffset(hitid) + layers * ((recHitTools_->zside(hitid) + 1) >> 1) - 1;
2137 std::unordered_map<DetId, const unsigned int>::const_iterator itcheck = hitMap.find(hitid);
2138 if (itcheck != hitMap.end()) {
2139 const HGCRecHit* hit = &(hits[itcheck->second]);
2140 sCEnergyOnLayer[scLayerId] += it_haf.second * hit->energy();
2141 }
2142 }
2143
2144 for (unsigned int layerId = 0; layerId < layers * 2; ++layerId) {
2145 if (!sCEnergyOnLayer[layerId])
2146 continue;
2147
2148 histograms.h_denom_simcluster_eta_perlayer[count].at(layerId)->Fill(sC[scId].eta());
2149 histograms.h_denom_simcluster_phi_perlayer[count].at(layerId)->Fill(sC[scId].phi());
2150
2151 if (lcsIt == lcsInSimClusterMap.end())
2152 continue;
2153 const auto& lcs = lcsIt->val;
2154
2155 auto getLCLayerId = [&](const unsigned int lcId) {
2156 const auto firstHitDetId = (clusters[lcId].hitsAndFractions())[0].first;
2157 if (recHitTools_->isBarrel(firstHitDetId))
2158 return 9999u;
2159 const unsigned int lcLayerId = recHitTools_->getLayerWithOffset(firstHitDetId) +
2160 layers * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
2161 return lcLayerId;
2162 };
2163
2164
2165 for (const auto& lcPair : lcs) {
2166 auto lcId = lcPair.first.index();
2167 if (mask[lcId] != 0.) {
2168 LogDebug("HGCalValidator") << "Skipping layer cluster " << lcId << " not belonging to mask" << std::endl;
2169 continue;
2170 }
2171
2172 if (getLCLayerId(lcId) != layerId)
2173 continue;
2174 histograms.h_score_simcluster2layercl_perlayer[count].at(layerId)->Fill(lcPair.second.second);
2175 histograms.h_sharedenergy_simcluster2layercl_perlayer[count].at(layerId)->Fill(
2176 lcPair.second.first / sCEnergyOnLayer[layerId], sCEnergyOnLayer[layerId]);
2177 histograms.h_energy_vs_score_simcluster2layercl_perlayer[count].at(layerId)->Fill(
2178 lcPair.second.second, lcPair.second.first / sCEnergyOnLayer[layerId]);
2179 }
2180 const auto assoc = std::count_if(std::begin(lcs), std::end(lcs), [&](const auto& obj) {
2181 if (getLCLayerId(obj.first.index()) != layerId)
2182 return false;
2183 else
2184 return obj.second.second < ScoreCutSCtoLC_;
2185 });
2186 if (assoc) {
2187 histograms.h_num_simcluster_eta_perlayer[count].at(layerId)->Fill(sC[scId].eta());
2188 histograms.h_num_simcluster_phi_perlayer[count].at(layerId)->Fill(sC[scId].phi());
2189 if (assoc > 1) {
2190 histograms.h_numDup_simcluster_eta_perlayer[count].at(layerId)->Fill(sC[scId].eta());
2191 histograms.h_numDup_simcluster_phi_perlayer[count].at(layerId)->Fill(sC[scId].phi());
2192 }
2193 const auto best = std::min_element(std::begin(lcs), std::end(lcs), [&](const auto& obj1, const auto& obj2) {
2194 if (getLCLayerId(obj1.first.index()) != layerId)
2195 return false;
2196 else if (getLCLayerId(obj2.first.index()) == layerId)
2197 return obj1.second.second < obj2.second.second;
2198 else
2199 return true;
2200 });
2201 histograms.h_sharedenergy_simcluster2layercl_vs_eta_perlayer[count].at(layerId)->Fill(
2202 sC[scId].eta(), best->second.first / sCEnergyOnLayer[layerId]);
2203 histograms.h_sharedenergy_simcluster2layercl_vs_phi_perlayer[count].at(layerId)->Fill(
2204 sC[scId].phi(), best->second.first / sCEnergyOnLayer[layerId]);
2205 }
2206 }
2207 }
2208 }
2209
2210 void HGVHistoProducerAlgo::fill_generic_cluster_histos(const Histograms& histograms,
2211 const int count,
2212 edm::Handle<reco::CaloClusterCollection> clusterHandle,
2213 const reco::CaloClusterCollection& clusters,
2214 edm::Handle<std::vector<CaloParticle>> caloParticleHandle,
2215 std::vector<CaloParticle> const& cP,
2216 std::vector<size_t> const& cPIndices,
2217 std::vector<size_t> const& cPSelectedIndices,
2218 std::unordered_map<DetId, const unsigned int> const& hitMap,
2219 std::map<double, double> cummatbudg,
2220 unsigned int layers,
2221 std::vector<int> thicknesses,
2222 const ticl::RecoToSimCollection& cpsInLayerClusterMap,
2223 const ticl::SimToRecoCollection& cPOnLayerMap,
2224 MultiVectorManager<HGCRecHit> const& hits) const {
2225
2226
2227
2228
2229
2230
2231
2232
2233 std::vector<int> tnlcpl(1000, 0);
2234
2235
2236 std::map<std::string, int> tnlcpthplus;
2237 tnlcpthplus.clear();
2238 std::map<std::string, int> tnlcpthminus;
2239 tnlcpthminus.clear();
2240
2241 for (std::vector<int>::iterator it = thicknesses.begin(); it != thicknesses.end(); ++it) {
2242 tnlcpthplus.insert(std::pair<std::string, int>(std::to_string(*it), 0));
2243 tnlcpthminus.insert(std::pair<std::string, int>(std::to_string(*it), 0));
2244 }
2245
2246 tnlcpthplus.insert(std::pair<std::string, int>("mixed", 0));
2247 tnlcpthminus.insert(std::pair<std::string, int>("mixed", 0));
2248
2249 layerClusters_to_CaloParticles(histograms,
2250 clusterHandle,
2251 clusters,
2252 caloParticleHandle,
2253 cP,
2254 cPIndices,
2255 cPSelectedIndices,
2256 hitMap,
2257 layers,
2258 cpsInLayerClusterMap,
2259 cPOnLayerMap,
2260 hits);
2261
2262
2263
2264 std::vector<double> tecpl(1000, 0.0);
2265
2266 std::vector<double> ldbar(1000, 0.0);
2267
2268
2269 double caloparteneplus = 0.;
2270 double caloparteneminus = 0.;
2271 for (const auto& cpId : cPIndices) {
2272 if (cP[cpId].eta() >= 0.) {
2273 caloparteneplus = caloparteneplus + cP[cpId].energy();
2274 } else if (cP[cpId].eta() < 0.) {
2275 caloparteneminus = caloparteneminus + cP[cpId].energy();
2276 }
2277 }
2278
2279
2280 for (const auto& lcId : clusters) {
2281 const auto seedid = lcId.seed();
2282 if (recHitTools_->isBarrel(seedid))
2283 continue;
2284 const double seedx = recHitTools_->getPosition(seedid).x();
2285 const double seedy = recHitTools_->getPosition(seedid).y();
2286 DetId maxid = findmaxhit(lcId, hitMap, hits);
2287
2288
2289 double maxx = recHitTools_->getPosition(maxid).x();
2290 double maxy = recHitTools_->getPosition(maxid).y();
2291
2292
2293 int nthhits120p = 0;
2294 int nthhits200p = 0;
2295 int nthhits300p = 0;
2296 int nthhitsscintp = 0;
2297 int nthhits120m = 0;
2298 int nthhits200m = 0;
2299 int nthhits300m = 0;
2300 int nthhitsscintm = 0;
2301
2302 double thickness = 0.;
2303
2304 int layerid = 0;
2305
2306
2307 int lay = 0;
2308
2309 std::string istr = "";
2310
2311 bool cluslay = true;
2312
2313 int zside = 0;
2314
2315 const auto& hits_and_fractions = lcId.hitsAndFractions();
2316 for (std::vector<std::pair<DetId, float>>::const_iterator it_haf = hits_and_fractions.begin();
2317 it_haf != hits_and_fractions.end();
2318 ++it_haf) {
2319 const DetId rh_detid = it_haf->first;
2320
2321 layerid = recHitTools_->getLayerWithOffset(rh_detid) + layers * ((recHitTools_->zside(rh_detid) + 1) >> 1) - 1;
2322 lay = recHitTools_->getLayerWithOffset(rh_detid);
2323 zside = recHitTools_->zside(rh_detid);
2324 if (rh_detid.det() == DetId::Forward || rh_detid.det() == DetId::HGCalEE || rh_detid.det() == DetId::HGCalHSi) {
2325 thickness = recHitTools_->getSiThickness(rh_detid);
2326 } else if (rh_detid.det() == DetId::HGCalHSc) {
2327 thickness = -1;
2328 } else {
2329 LogDebug("HGCalValidator") << "These are HGCal layer clusters, you shouldn't be here !!! " << layerid << "\n";
2330 continue;
2331 }
2332
2333
2334 std::string curistr = std::to_string((int)thickness);
2335 std::string lay_string = std::to_string(layerid);
2336 while (lay_string.size() < 2)
2337 lay_string.insert(0, "0");
2338 curistr += "_" + lay_string;
2339 if (cluslay) {
2340 tnlcpl[layerid]++;
2341 istr = curistr;
2342 cluslay = false;
2343 }
2344
2345 if ((thickness == 120.) && (recHitTools_->zside(rh_detid) > 0.)) {
2346 nthhits120p++;
2347 } else if ((thickness == 120.) && (recHitTools_->zside(rh_detid) < 0.)) {
2348 nthhits120m++;
2349 } else if ((thickness == 200.) && (recHitTools_->zside(rh_detid) > 0.)) {
2350 nthhits200p++;
2351 } else if ((thickness == 200.) && (recHitTools_->zside(rh_detid) < 0.)) {
2352 nthhits200m++;
2353 } else if ((thickness == 300.) && (recHitTools_->zside(rh_detid) > 0.)) {
2354 nthhits300p++;
2355 } else if ((thickness == 300.) && (recHitTools_->zside(rh_detid) < 0.)) {
2356 nthhits300m++;
2357 } else if ((thickness == -1) && (recHitTools_->zside(rh_detid) > 0.)) {
2358 nthhitsscintp++;
2359 } else if ((thickness == -1) && (recHitTools_->zside(rh_detid) < 0.)) {
2360 nthhitsscintm++;
2361 } else {
2362 LogDebug("HGCalValidator")
2363 << " You are running a geometry that contains thicknesses different than the normal ones. "
2364 << "\n";
2365 }
2366
2367 std::unordered_map<DetId, const unsigned int>::const_iterator itcheck = hitMap.find(rh_detid);
2368 if (itcheck == hitMap.end()) {
2369 std::ostringstream st1;
2370 if ((rh_detid.det() == DetId::HGCalEE) || (rh_detid.det() == DetId::HGCalHSi)) {
2371 st1 << HGCSiliconDetId(rh_detid);
2372 } else if (rh_detid.det() == DetId::HGCalHSc) {
2373 st1 << HGCScintillatorDetId(rh_detid);
2374 } else {
2375 st1 << HFNoseDetId(rh_detid);
2376 }
2377 LogDebug("HGCalValidator") << " You shouldn't be here - Unable to find a hit " << rh_detid.rawId() << " "
2378 << rh_detid.det() << " " << st1.str() << "\n";
2379 continue;
2380 }
2381
2382 const HGCRecHit* hit = &(hits[itcheck->second]);
2383
2384
2385 const double hit_x = recHitTools_->getPosition(rh_detid).x();
2386 const double hit_y = recHitTools_->getPosition(rh_detid).y();
2387 double distancetoseed = distance(seedx, seedy, hit_x, hit_y);
2388 double distancetomax = distance(maxx, maxy, hit_x, hit_y);
2389 if (distancetoseed != 0. && histograms.h_distancetoseedcell_perthickperlayer.count(curistr)) {
2390 histograms.h_distancetoseedcell_perthickperlayer.at(curistr)->Fill(distancetoseed);
2391 }
2392
2393 if (distancetoseed != 0. && histograms.h_distancetoseedcell_perthickperlayer_eneweighted.count(curistr)) {
2394 histograms.h_distancetoseedcell_perthickperlayer_eneweighted.at(curistr)->Fill(distancetoseed, hit->energy());
2395 }
2396
2397 if (distancetomax != 0. && histograms.h_distancetomaxcell_perthickperlayer.count(curistr)) {
2398 histograms.h_distancetomaxcell_perthickperlayer.at(curistr)->Fill(distancetomax);
2399 }
2400
2401 if (distancetomax != 0. && histograms.h_distancetomaxcell_perthickperlayer_eneweighted.count(curistr)) {
2402 histograms.h_distancetomaxcell_perthickperlayer_eneweighted.at(curistr)->Fill(distancetomax, hit->energy());
2403 }
2404
2405 }
2406
2407
2408 if ((nthhits120p != 0 && nthhits200p != 0) || (nthhits120p != 0 && nthhits300p != 0) ||
2409 (nthhits120p != 0 && nthhitsscintp != 0) || (nthhits200p != 0 && nthhits300p != 0) ||
2410 (nthhits200p != 0 && nthhitsscintp != 0) || (nthhits300p != 0 && nthhitsscintp != 0)) {
2411 tnlcpthplus["mixed"]++;
2412 } else if ((nthhits120p != 0 || nthhits200p != 0 || nthhits300p != 0 || nthhitsscintp != 0)) {
2413
2414 tnlcpthplus[std::to_string((int)thickness)]++;
2415 }
2416 if ((nthhits120m != 0 && nthhits200m != 0) || (nthhits120m != 0 && nthhits300m != 0) ||
2417 (nthhits120m != 0 && nthhitsscintm != 0) || (nthhits200m != 0 && nthhits300m != 0) ||
2418 (nthhits200m != 0 && nthhitsscintm != 0) || (nthhits300m != 0 && nthhitsscintm != 0)) {
2419 tnlcpthminus["mixed"]++;
2420 } else if ((nthhits120m != 0 || nthhits200m != 0 || nthhits300m != 0 || nthhitsscintm != 0)) {
2421
2422 tnlcpthminus[std::to_string((int)thickness)]++;
2423 }
2424
2425
2426 std::vector<int> bigamoth;
2427 bigamoth.clear();
2428 if (zside > 0) {
2429 bigamoth.push_back(nthhits120p);
2430 bigamoth.push_back(nthhits200p);
2431 bigamoth.push_back(nthhits300p);
2432 bigamoth.push_back(nthhitsscintp);
2433 } else if (zside < 0) {
2434 bigamoth.push_back(nthhits120m);
2435 bigamoth.push_back(nthhits200m);
2436 bigamoth.push_back(nthhits300m);
2437 bigamoth.push_back(nthhitsscintm);
2438 }
2439 auto bgth = std::max_element(bigamoth.begin(), bigamoth.end());
2440 istr = std::to_string(thicknesses[std::distance(bigamoth.begin(), bgth)]);
2441 std::string lay_string = std::to_string(layerid);
2442 while (lay_string.size() < 2)
2443 lay_string.insert(0, "0");
2444 istr += "_" + lay_string;
2445
2446
2447 if (histograms.h_cellsnum_perthickperlayer.count(istr)) {
2448 histograms.h_cellsnum_perthickperlayer.at(istr)->Fill(hits_and_fractions.size());
2449 }
2450
2451
2452 double distancebetseedandmax = distance(seedx, seedy, maxx, maxy);
2453
2454 std::string seedstr = std::to_string((int)recHitTools_->getSiThickness(seedid)) + "_" + std::to_string(layerid);
2455 seedstr += "_" + lay_string;
2456 if (histograms.h_distancebetseedandmaxcell_perthickperlayer.count(seedstr)) {
2457 histograms.h_distancebetseedandmaxcell_perthickperlayer.at(seedstr)->Fill(distancebetseedandmax);
2458 }
2459 const auto lc_en = lcId.energy();
2460 if (histograms.h_distancebetseedandmaxcellvsclusterenergy_perthickperlayer.count(seedstr)) {
2461 histograms.h_distancebetseedandmaxcellvsclusterenergy_perthickperlayer.at(seedstr)->Fill(distancebetseedandmax,
2462 lc_en);
2463 }
2464
2465
2466 tecpl[layerid] = tecpl[layerid] + lc_en;
2467 ldbar[layerid] = ldbar[layerid] + lc_en * cummatbudg[(double)lay];
2468
2469 }
2470
2471
2472 double sumeneallcluspl = 0.;
2473 double sumeneallclusmi = 0.;
2474
2475 double sumldbarpl = 0.;
2476 double sumldbarmi = 0.;
2477
2478 for (unsigned ilayer = 0; ilayer < layers * 2; ++ilayer) {
2479 if (histograms.h_clusternum_perlayer.count(ilayer)) {
2480 histograms.h_clusternum_perlayer.at(ilayer)->Fill(tnlcpl[ilayer]);
2481 }
2482
2483
2484 if (ilayer < layers) {
2485 if (histograms.h_energyclustered_perlayer.count(ilayer)) {
2486 if (caloparteneminus != 0.) {
2487 histograms.h_energyclustered_perlayer.at(ilayer)->Fill(100. * tecpl[ilayer] / caloparteneminus);
2488 }
2489 }
2490
2491 sumeneallclusmi = sumeneallclusmi + tecpl[ilayer];
2492
2493 sumldbarmi = sumldbarmi + ldbar[ilayer];
2494 } else {
2495 if (histograms.h_energyclustered_perlayer.count(ilayer)) {
2496 if (caloparteneplus != 0.) {
2497 histograms.h_energyclustered_perlayer.at(ilayer)->Fill(100. * tecpl[ilayer] / caloparteneplus);
2498 }
2499 }
2500
2501 sumeneallcluspl = sumeneallcluspl + tecpl[ilayer];
2502
2503 sumldbarpl = sumldbarpl + ldbar[ilayer];
2504 }
2505
2506 }
2507
2508
2509 for (std::vector<int>::iterator it = thicknesses.begin(); it != thicknesses.end(); ++it) {
2510 if (histograms.h_clusternum_perthick.count(*it)) {
2511 histograms.h_clusternum_perthick.at(*it)->Fill(tnlcpthplus[std::to_string(*it)]);
2512 histograms.h_clusternum_perthick.at(*it)->Fill(tnlcpthminus[std::to_string(*it)]);
2513 }
2514 }
2515
2516 histograms.h_mixedhitscluster_zplus[count]->Fill(tnlcpthplus["mixed"]);
2517 histograms.h_mixedhitscluster_zminus[count]->Fill(tnlcpthminus["mixed"]);
2518
2519
2520 if (caloparteneplus != 0.) {
2521 histograms.h_energyclustered_zplus[count]->Fill(100. * sumeneallcluspl / caloparteneplus);
2522 }
2523 if (caloparteneminus != 0.) {
2524 histograms.h_energyclustered_zminus[count]->Fill(100. * sumeneallclusmi / caloparteneminus);
2525 }
2526
2527
2528 histograms.h_longdepthbarycentre_zplus[count]->Fill(sumldbarpl / sumeneallcluspl);
2529 histograms.h_longdepthbarycentre_zminus[count]->Fill(sumldbarmi / sumeneallclusmi);
2530 }
2531
2532 void HGVHistoProducerAlgo::tracksters_to_SimTracksters_fp(const Histograms& histograms,
2533 const int count,
2534 const TracksterToTracksterMap& trackstersToSimTrackstersMap,
2535 const TracksterToTracksterMap& simTrackstersToTrackstersMap,
2536 const validationType valType,
2537 const SimClusterToCaloParticleMap& scToCpMap,
2538 const std::vector<size_t>& cPIndices,
2539 const std::vector<size_t>& cPSelectedIndices,
2540 const edm::ProductID& cPHandle_id) const {
2541 const auto nTracksters = trackstersToSimTrackstersMap.getMap().size();
2542 const auto nSimTracksters = simTrackstersToTrackstersMap.getMap().size();
2543 std::vector<int> tracksters_FakeMerge(nTracksters, 0);
2544 std::vector<int> tracksters_PurityDuplicate(nSimTracksters, 0);
2545 auto getCPId = [](const ticl::Trackster& simTS,
2546 const edm::ProductID& cPHandle_id,
2547 const SimClusterToCaloParticleMap& scToCpMap) {
2548 const auto productID = simTS.seedID();
2549 if (productID == cPHandle_id) {
2550 return simTS.seedIndex();
2551 } else {
2552 return int(scToCpMap.at(simTS.seedIndex()).index());
2553 }
2554 };
2555
2556 auto ScoreCutSTStoTSPurDup = ScoreCutSTStoTSPurDup_[0];
2557 auto ScoreCutTStoSTSFakeMerge = ScoreCutTStoSTSFakeMerge_[0];
2558 for (unsigned int tracksterIndex = 0; tracksterIndex < nTracksters; ++tracksterIndex) {
2559 const auto& trackster = *(trackstersToSimTrackstersMap.getRefFirst(tracksterIndex));
2560 if (trackster.vertices().empty())
2561 continue;
2562 float iTS_eta = trackster.barycenter().eta();
2563 float iTS_phi = trackster.barycenter().phi();
2564 float iTS_en = trackster.raw_energy();
2565 float iTS_pt = trackster.raw_pt();
2566 histograms.h_denom_trackster_eta[valType][count]->Fill(iTS_eta);
2567 histograms.h_denom_trackster_phi[valType][count]->Fill(iTS_phi);
2568 histograms.h_denom_trackster_en[valType][count]->Fill(iTS_en);
2569 histograms.h_denom_trackster_pt[valType][count]->Fill(iTS_pt);
2570
2571
2572 for (unsigned int i = 0; i < trackstersToSimTrackstersMap[tracksterIndex].size(); ++i) {
2573 auto sharedEnergy = trackstersToSimTrackstersMap[tracksterIndex][i].sharedEnergy();
2574 auto score = trackstersToSimTrackstersMap[tracksterIndex][i].score();
2575 float sharedEnergyFraction = sharedEnergy / trackster.raw_energy();
2576 if (i == 0) {
2577 histograms.h_score_trackster2bestCaloparticle[valType][count]->Fill(score);
2578 histograms.h_sharedenergy_trackster2bestCaloparticle[valType][count]->Fill(sharedEnergyFraction);
2579 histograms.h_sharedenergy_trackster2bestCaloparticle_vs_eta[valType][count]->Fill(trackster.barycenter().eta(),
2580 sharedEnergy);
2581 histograms.h_sharedenergy_trackster2bestCaloparticle_vs_phi[valType][count]->Fill(trackster.barycenter().phi(),
2582 sharedEnergy);
2583 histograms.h_energy_vs_score_trackster2bestCaloparticle[valType][count]->Fill(score, sharedEnergyFraction);
2584 }
2585 if (i == 1) {
2586 histograms.h_score_trackster2bestCaloparticle2[valType][count]->Fill(score);
2587 histograms.h_sharedenergy_trackster2bestCaloparticle2[valType][count]->Fill(sharedEnergyFraction);
2588 histograms.h_energy_vs_score_trackster2bestCaloparticle2[valType][count]->Fill(score, sharedEnergyFraction);
2589 }
2590 histograms.h_score_trackster2caloparticle[valType][count]->Fill(score);
2591 histograms.h_sharedenergy_trackster2caloparticle[valType][count]->Fill(sharedEnergyFraction);
2592 histograms.h_energy_vs_score_trackster2caloparticle[valType][count]->Fill(score, sharedEnergyFraction);
2593 tracksters_FakeMerge[tracksterIndex] += score < ScoreCutTStoSTSFakeMerge;
2594 }
2595
2596 if (tracksters_FakeMerge[tracksterIndex] > 0) {
2597 histograms.h_num_trackster_eta[valType][count]->Fill(iTS_eta);
2598 histograms.h_num_trackster_phi[valType][count]->Fill(iTS_phi);
2599 histograms.h_num_trackster_en[valType][count]->Fill(iTS_en);
2600 histograms.h_num_trackster_pt[valType][count]->Fill(iTS_pt);
2601
2602 if (tracksters_FakeMerge[tracksterIndex] > 1) {
2603 histograms.h_numMerge_trackster_eta[valType][count]->Fill(iTS_eta);
2604 histograms.h_numMerge_trackster_phi[valType][count]->Fill(iTS_phi);
2605 histograms.h_numMerge_trackster_en[valType][count]->Fill(iTS_en);
2606 histograms.h_numMerge_trackster_pt[valType][count]->Fill(iTS_pt);
2607 }
2608 }
2609 }
2610
2611
2612
2613
2614 for (unsigned int simTracksterIndex = 0; simTracksterIndex < nSimTracksters; ++simTracksterIndex) {
2615 const auto& simTrackster = *(simTrackstersToTrackstersMap.getRefFirst(simTracksterIndex));
2616 const auto cpId = getCPId(simTrackster, cPHandle_id, scToCpMap);
2617 if (std::find(cPSelectedIndices.begin(), cPSelectedIndices.end(), cpId) == cPSelectedIndices.end())
2618 continue;
2619 const auto sts_eta = simTrackster.barycenter().eta();
2620 const auto sts_phi = simTrackster.barycenter().phi();
2621 const auto sts_en = simTrackster.raw_energy();
2622 const auto sts_pt = simTrackster.raw_pt();
2623 float inv_simtrackster_energy = 1.f / sts_en;
2624 histograms.h_denom_caloparticle_eta[valType][count]->Fill(sts_eta);
2625 histograms.h_denom_caloparticle_phi[valType][count]->Fill(sts_phi);
2626 histograms.h_denom_caloparticle_en[valType][count]->Fill(sts_en);
2627 histograms.h_denom_caloparticle_pt[valType][count]->Fill(sts_pt);
2628
2629
2630
2631
2632
2633
2634 bool sts_considered_efficient = false;
2635 bool sts_considered_pure = false;
2636
2637 for (unsigned int i = 0; i < simTrackstersToTrackstersMap[simTracksterIndex].size(); ++i) {
2638 const auto sharedEnergy = simTrackstersToTrackstersMap[simTracksterIndex][i].sharedEnergy();
2639 const auto score = simTrackstersToTrackstersMap[simTracksterIndex][i].score();
2640 float sharedEnergyFraction = sharedEnergy * inv_simtrackster_energy;
2641 if (i == 0) {
2642 histograms.h_scorePur_caloparticle2trackster[valType][count]->Fill(score);
2643 histograms.h_sharedenergy_caloparticle2trackster_assoc[valType][count]->Fill(sharedEnergyFraction);
2644 histograms.h_energy_vs_score_caloparticle2bestTrackster[valType][count]->Fill(score, sharedEnergyFraction);
2645 histograms.h_sharedenergy_caloparticle2trackster_assoc_vs_eta[valType][count]->Fill(sts_eta,
2646 sharedEnergyFraction);
2647 histograms.h_sharedenergy_caloparticle2trackster_assoc_vs_phi[valType][count]->Fill(sts_phi,
2648 sharedEnergyFraction);
2649 }
2650
2651 if (i == 1) {
2652 histograms.h_scoreDupl_caloparticle2trackster[valType][count]->Fill(score);
2653 histograms.h_sharedenergy_caloparticle2trackster_assoc2[valType][count]->Fill(sharedEnergyFraction);
2654 histograms.h_energy_vs_score_caloparticle2bestTrackster2[valType][count]->Fill(score, sharedEnergyFraction);
2655 }
2656
2657 histograms.h_score_caloparticle2trackster[valType][count]->Fill(score);
2658 histograms.h_sharedenergy_caloparticle2trackster[valType][count]->Fill(sharedEnergyFraction);
2659 histograms.h_energy_vs_score_caloparticle2trackster[valType][count]->Fill(score, sharedEnergyFraction);
2660
2661
2662 if (!sts_considered_efficient && (sharedEnergyFraction >= minTSTSharedEneFracEfficiency_)) {
2663 sts_considered_efficient = true;
2664 histograms.h_numEff_caloparticle_eta[valType][count]->Fill(sts_eta);
2665 histograms.h_numEff_caloparticle_phi[valType][count]->Fill(sts_phi);
2666 histograms.h_numEff_caloparticle_en[valType][count]->Fill(sts_en);
2667 histograms.h_numEff_caloparticle_pt[valType][count]->Fill(sts_pt);
2668 }
2669
2670 if (score < ScoreCutSTStoTSPurDup) {
2671 if (tracksters_PurityDuplicate[simTracksterIndex] < 1)
2672 tracksters_PurityDuplicate[simTracksterIndex]++;
2673 if (sts_considered_pure)
2674 tracksters_PurityDuplicate[simTracksterIndex]++;
2675 sts_considered_pure = true;
2676 }
2677
2678 }
2679 if (tracksters_PurityDuplicate[simTracksterIndex] > 0) {
2680 histograms.h_num_caloparticle_eta[valType][count]->Fill(sts_eta);
2681 histograms.h_num_caloparticle_phi[valType][count]->Fill(sts_phi);
2682 histograms.h_num_caloparticle_en[valType][count]->Fill(sts_en);
2683 histograms.h_num_caloparticle_pt[valType][count]->Fill(sts_pt);
2684
2685 if (tracksters_PurityDuplicate[simTracksterIndex] > 1) {
2686 histograms.h_numDup_trackster_eta[valType][count]->Fill(sts_eta);
2687 histograms.h_numDup_trackster_phi[valType][count]->Fill(sts_phi);
2688 histograms.h_numDup_trackster_en[valType][count]->Fill(sts_en);
2689 histograms.h_numDup_trackster_pt[valType][count]->Fill(sts_pt);
2690 }
2691 }
2692
2693 }
2694 }
2695
2696 void HGVHistoProducerAlgo::fill_trackster_histos(
2697 const Histograms& histograms,
2698 const int count,
2699 const ticl::TracksterCollection& tracksters,
2700 const reco::CaloClusterCollection& layerClusters,
2701 const ticl::TracksterCollection& simTSs,
2702 const ticl::TracksterCollection& simTSs_fromCP,
2703 const std::map<unsigned int, std::vector<unsigned int>>& cpToSc_SimTrackstersMap,
2704 std::vector<SimCluster> const& sC,
2705 const edm::ProductID& cPHandle_id,
2706 std::vector<CaloParticle> const& cP,
2707 std::vector<size_t> const& cPIndices,
2708 std::vector<size_t> const& cPSelectedIndices,
2709 std::unordered_map<DetId, const unsigned int> const& hitMap,
2710 unsigned int layers,
2711 MultiVectorManager<HGCRecHit> const& hits,
2712 bool mapsFound,
2713 const edm::Handle<TracksterToTracksterMap>& trackstersToSimTrackstersByLCsMapH,
2714 const edm::Handle<TracksterToTracksterMap>& simTrackstersToTrackstersByLCsMapH,
2715 const edm::Handle<TracksterToTracksterMap>& trackstersToSimTrackstersFromCPsByLCsMapH,
2716 const edm::Handle<TracksterToTracksterMap>& simTrackstersFromCPsToTrackstersByLCsMapH,
2717 const edm::Handle<TracksterToTracksterMap>& trackstersToSimTrackstersByHitsMapH,
2718 const edm::Handle<TracksterToTracksterMap>& simTrackstersToTrackstersByHitsMapH,
2719 const edm::Handle<TracksterToTracksterMap>& trackstersToSimTrackstersFromCPsByHitsMapH,
2720 const edm::Handle<TracksterToTracksterMap>& simTrackstersFromCPsToTrackstersByHitsMapH,
2721 const SimClusterToCaloParticleMap& scToCpMap) const {
2722
2723
2724
2725
2726 int totNTstZm = 0;
2727 int totNTstZp = 0;
2728
2729 int totNContTstZp = 0;
2730 int totNContTstZm = 0;
2731
2732 int totNNotContTstZp = 0;
2733 int totNNotContTstZm = 0;
2734
2735 std::vector<bool> contTracksters;
2736 contTracksters.clear();
2737
2738
2739 std::unordered_map<unsigned int, std::vector<unsigned int>> multiplicity;
2740
2741 std::unordered_map<unsigned int, std::vector<unsigned int>> multiplicity_vs_layer;
2742
2743
2744
2745
2746
2747
2748 const auto nTracksters = tracksters.size();
2749
2750 for (unsigned int tstId = 0; tstId < nTracksters; ++tstId) {
2751 const auto& tst = tracksters[tstId];
2752 if (tst.vertices().empty())
2753 continue;
2754
2755 if (tst.barycenter().z() < 0.)
2756 totNTstZm++;
2757 else if (tst.barycenter().z() > 0.)
2758 totNTstZp++;
2759
2760
2761 int tnLcInTst = 0;
2762
2763
2764
2765 std::vector<int> tnLcInTstperlay(1000, 0);
2766
2767
2768
2769 std::set<unsigned int> trackster_layers;
2770
2771 bool tracksterInZplus = false;
2772 bool tracksterInZminus = false;
2773
2774
2775 for (const auto lcId : tst.vertices()) {
2776
2777 const auto& hits_and_fractions = layerClusters[lcId].hitsAndFractions();
2778 if (recHitTools_->isBarrel(hits_and_fractions[0].first))
2779 continue;
2780
2781 multiplicity[tstId].emplace_back(hits_and_fractions.size());
2782
2783 const auto firstHitDetId = hits_and_fractions[0].first;
2784 if (recHitTools_->isBarrel(firstHitDetId))
2785 continue;
2786
2787 const auto layerid = recHitTools_->getLayerWithOffset(firstHitDetId) +
2788 layers * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
2789 trackster_layers.insert(layerid);
2790 multiplicity_vs_layer[tstId].emplace_back(layerid);
2791
2792 tnLcInTstperlay[layerid]++;
2793 tnLcInTst++;
2794
2795 if (recHitTools_->zside(firstHitDetId) > 0.)
2796 tracksterInZplus = true;
2797 else if (recHitTools_->zside(firstHitDetId) < 0.)
2798 tracksterInZminus = true;
2799 }
2800
2801
2802 for (unsigned ilayer = 0; ilayer < layers * 2; ++ilayer) {
2803 if (histograms.h_clusternum_in_trackster_perlayer[count].count(ilayer) && tnLcInTstperlay[ilayer] != 0) {
2804 histograms.h_clusternum_in_trackster_perlayer[count].at(ilayer)->Fill((float)tnLcInTstperlay[ilayer]);
2805 }
2806
2807 if (tnLcInTstperlay[ilayer] != 0) {
2808 histograms.h_clusternum_in_trackster_vs_layer[count]->Fill((float)ilayer, (float)tnLcInTstperlay[ilayer]);
2809 }
2810 }
2811
2812
2813 std::vector<int> trackster_layers_vec(trackster_layers.begin(), trackster_layers.end());
2814
2815 bool contiTrackster = false;
2816
2817 if (trackster_layers_vec.size() >= 3) {
2818 for (unsigned int iLayer = 1; iLayer < trackster_layers_vec.size() - 1; ++iLayer) {
2819 if ((trackster_layers_vec[iLayer - 1] + 1 == trackster_layers_vec[iLayer]) &&
2820 (trackster_layers_vec[iLayer + 1] - 1 == trackster_layers_vec[iLayer])) {
2821
2822 if (tracksterInZplus)
2823 totNContTstZp++;
2824 else if (tracksterInZminus)
2825 totNContTstZm++;
2826
2827 contiTrackster = true;
2828 break;
2829 }
2830 }
2831 }
2832
2833 if (!contiTrackster) {
2834 if (tracksterInZplus)
2835 totNNotContTstZp++;
2836 else if (tracksterInZminus)
2837 totNNotContTstZm++;
2838 }
2839
2840
2841 contTracksters.push_back(contiTrackster);
2842
2843 histograms.h_clusternum_in_trackster[count]->Fill(tnLcInTst);
2844
2845 for (unsigned int lc = 0; lc < multiplicity[tstId].size(); ++lc) {
2846
2847 float mlp = std::count(std::begin(multiplicity[tstId]), std::end(multiplicity[tstId]), multiplicity[tstId][lc]);
2848
2849
2850 histograms.h_multiplicityOfLCinTST[count]->Fill(mlp, multiplicity[tstId][lc]);
2851
2852
2853 histograms.h_multiplicity_numberOfEventsHistogram[count]->Fill(mlp);
2854
2855
2856 if (multiplicity_vs_layer[tstId][lc] < layers) {
2857 histograms.h_multiplicityOfLCinTST_vs_layercluster_zminus[count]->Fill(mlp, multiplicity_vs_layer[tstId][lc]);
2858 histograms.h_multiplicity_zminus_numberOfEventsHistogram[count]->Fill(mlp);
2859 } else {
2860 histograms.h_multiplicityOfLCinTST_vs_layercluster_zplus[count]->Fill(
2861 mlp, multiplicity_vs_layer[tstId][lc] - layers);
2862 histograms.h_multiplicity_zplus_numberOfEventsHistogram[count]->Fill(mlp);
2863 }
2864
2865 histograms.h_multiplicityOfLCinTST_vs_layerclusterenergy[count]->Fill(mlp,
2866 layerClusters[tst.vertices(lc)].energy());
2867 }
2868
2869 if (!trackster_layers.empty()) {
2870 histograms.h_trackster_x[count]->Fill(tst.barycenter().x());
2871 histograms.h_trackster_y[count]->Fill(tst.barycenter().y());
2872 histograms.h_trackster_z[count]->Fill(tst.barycenter().z());
2873 histograms.h_trackster_eta[count]->Fill(tst.barycenter().eta());
2874 histograms.h_trackster_phi[count]->Fill(tst.barycenter().phi());
2875
2876 histograms.h_trackster_firstlayer[count]->Fill((float)*trackster_layers.begin());
2877 histograms.h_trackster_lastlayer[count]->Fill((float)*trackster_layers.rbegin());
2878 histograms.h_trackster_layersnum[count]->Fill((float)trackster_layers.size());
2879
2880 histograms.h_trackster_pt[count]->Fill(tst.raw_pt());
2881 histograms.h_trackster_energy[count]->Fill(tst.raw_energy());
2882 }
2883
2884 }
2885
2886 histograms.h_tracksternum[count]->Fill(totNTstZm + totNTstZp);
2887 histograms.h_conttracksternum[count]->Fill(totNContTstZp + totNContTstZm);
2888 histograms.h_nonconttracksternum[count]->Fill(totNNotContTstZp + totNNotContTstZm);
2889 if (mapsFound) {
2890 const auto& trackstersToSimTrackstersByLCsMap = *trackstersToSimTrackstersByLCsMapH;
2891 const auto& simTrackstersToTrackstersByLCsMap = *simTrackstersToTrackstersByLCsMapH;
2892 const auto& trackstersToSimTrackstersFromCPsByLCsMap = *trackstersToSimTrackstersFromCPsByLCsMapH;
2893 const auto& simTrackstersFromCPsToTrackstersByLCsMap = *simTrackstersFromCPsToTrackstersByLCsMapH;
2894 const auto& trackstersToSimTrackstersByHitsMap = *trackstersToSimTrackstersByHitsMapH;
2895 const auto& simTrackstersToTrackstersByHitsMap = *simTrackstersToTrackstersByHitsMapH;
2896 const auto& trackstersToSimTrackstersFromCPsByHitsMap = *trackstersToSimTrackstersFromCPsByHitsMapH;
2897 const auto& simTrackstersFromCPsToTrackstersByHitsMap = *simTrackstersFromCPsToTrackstersByHitsMapH;
2898
2899 tracksters_to_SimTracksters_fp(histograms,
2900 count,
2901 trackstersToSimTrackstersByLCsMap,
2902 simTrackstersToTrackstersByLCsMap,
2903 validationType::byLCs,
2904 scToCpMap,
2905 cPIndices,
2906 cPSelectedIndices,
2907 cPHandle_id);
2908
2909 tracksters_to_SimTracksters_fp(histograms,
2910 count,
2911 trackstersToSimTrackstersFromCPsByLCsMap,
2912 simTrackstersFromCPsToTrackstersByLCsMap,
2913 validationType::byLCs_CP,
2914 scToCpMap,
2915 cPIndices,
2916 cPSelectedIndices,
2917 cPHandle_id);
2918
2919 tracksters_to_SimTracksters_fp(histograms,
2920 count,
2921 trackstersToSimTrackstersFromCPsByHitsMap,
2922 simTrackstersFromCPsToTrackstersByHitsMap,
2923 validationType::byHits_CP,
2924 scToCpMap,
2925 cPIndices,
2926 cPSelectedIndices,
2927 cPHandle_id);
2928
2929 tracksters_to_SimTracksters_fp(histograms,
2930 count,
2931 trackstersToSimTrackstersByHitsMap,
2932 simTrackstersToTrackstersByHitsMap,
2933 validationType::byHits,
2934 scToCpMap,
2935 cPIndices,
2936 cPSelectedIndices,
2937 cPHandle_id);
2938 }
2939 }
2940
2941 double HGVHistoProducerAlgo::distance2(const double x1,
2942 const double y1,
2943 const double x2,
2944 const double y2) const {
2945 const double dx = x1 - x2;
2946 const double dy = y1 - y2;
2947 return (dx * dx + dy * dy);
2948 }
2949 double HGVHistoProducerAlgo::distance(const double x1,
2950 const double y1,
2951 const double x2,
2952 const double y2) const {
2953 return std::sqrt(distance2(x1, y1, x2, y2));
2954 }
2955
2956 void HGVHistoProducerAlgo::setRecHitTools(std::shared_ptr<hgcal::RecHitTools> recHitTools) {
2957 recHitTools_ = recHitTools;
2958 }
2959
2960 DetId HGVHistoProducerAlgo::findmaxhit(const reco::CaloCluster& cluster,
2961 std::unordered_map<DetId, const unsigned int> const& hitMap,
2962 MultiVectorManager<HGCRecHit> const& hits) const {
2963 const auto& hits_and_fractions = cluster.hitsAndFractions();
2964
2965 DetId themaxid;
2966 double maxene = 0.;
2967 for (std::vector<std::pair<DetId, float>>::const_iterator it_haf = hits_and_fractions.begin();
2968 it_haf != hits_and_fractions.end();
2969 ++it_haf) {
2970 const DetId rh_detid = it_haf->first;
2971 if (recHitTools_->isBarrel(rh_detid))
2972 continue;
2973 const auto hitEn = (hits[hitMap.find(rh_detid)->second]).energy();
2974 if (maxene < hitEn) {
2975 maxene = hitEn;
2976 themaxid = rh_detid;
2977 }
2978 }
2979
2980 return themaxid;
2981 }
2982
2983 double HGVHistoProducerAlgo::getEta(double eta) const {
2984 if (useFabsEta_)
2985 return fabs(eta);
2986 else
2987 return eta;
2988 }
