CMS 3D CMS Logo
HGVHistoProducerAlgo.cc
Go to the documentation of this file.
1 #include <numeric>
2 #include <iomanip>
3 #include <sstream>
4 
5 #include "Validation/HGCalValidation/interface/HGVHistoProducerAlgo.h"
6 #include "FWCore/MessageLogger/interface/MessageLogger.h"
7 #include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
8 #include "TMath.h"
9 #include "TLatex.h"
10 #include "TF1.h"
11 
12 using namespace std;
13 
14 //Parameters for the score cut. Later, this will become part of the
15 //configuration parameter for the HGCAL associator.
16 const double ScoreCutLCtoCP_ = 0.1;
17 const double ScoreCutCPtoLC_ = 0.1;
18 const double ScoreCutLCtoSC_ = 0.1;
19 const double ScoreCutSCtoLC_ = 0.1;
20 const double ScoreCutTStoSTSFakeMerge_[] = {0.6, FLT_MIN}; //1.e-09
21 const double ScoreCutSTStoTSPurDup_[] = {0.2, FLT_MIN}; //1.e-11
22 
23 HGVHistoProducerAlgo::HGVHistoProducerAlgo(const edm::ParameterSet& pset)
24  : //parameters for eta
25  minEta_(pset.getParameter<double>("minEta")),
26  maxEta_(pset.getParameter<double>("maxEta")),
27  nintEta_(pset.getParameter<int>("nintEta")),
28  useFabsEta_(pset.getParameter<bool>("useFabsEta")),
29 
30  //parameters for energy
31  minEne_(pset.getParameter<double>("minEne")),
32  maxEne_(pset.getParameter<double>("maxEne")),
33  nintEne_(pset.getParameter<int>("nintEne")),
34 
35  //parameters for pt
36  minPt_(pset.getParameter<double>("minPt")),
37  maxPt_(pset.getParameter<double>("maxPt")),
38  nintPt_(pset.getParameter<int>("nintPt")),
39 
40  //parameters for phi
41  minPhi_(pset.getParameter<double>("minPhi")),
42  maxPhi_(pset.getParameter<double>("maxPhi")),
43  nintPhi_(pset.getParameter<int>("nintPhi")),
44 
45  //parameters for counting mixed hits SimClusters
46  minMixedHitsSimCluster_(pset.getParameter<double>("minMixedHitsSimCluster")),
47  maxMixedHitsSimCluster_(pset.getParameter<double>("maxMixedHitsSimCluster")),
48  nintMixedHitsSimCluster_(pset.getParameter<int>("nintMixedHitsSimCluster")),
49 
50  //parameters for counting mixed hits clusters
51  minMixedHitsCluster_(pset.getParameter<double>("minMixedHitsCluster")),
52  maxMixedHitsCluster_(pset.getParameter<double>("maxMixedHitsCluster")),
53  nintMixedHitsCluster_(pset.getParameter<int>("nintMixedHitsCluster")),
54 
55  //parameters for the total amount of energy clustered by all layer clusters (fraction over CaloParticless)
56  minEneCl_(pset.getParameter<double>("minEneCl")),
57  maxEneCl_(pset.getParameter<double>("maxEneCl")),
58  nintEneCl_(pset.getParameter<int>("nintEneCl")),
59 
60  //parameters for the longitudinal depth barycenter.
61  minLongDepBary_(pset.getParameter<double>("minLongDepBary")),
62  maxLongDepBary_(pset.getParameter<double>("maxLongDepBary")),
63  nintLongDepBary_(pset.getParameter<int>("nintLongDepBary")),
64 
65  //parameters for z positionof vertex plots
66  minZpos_(pset.getParameter<double>("minZpos")),
67  maxZpos_(pset.getParameter<double>("maxZpos")),
68  nintZpos_(pset.getParameter<int>("nintZpos")),
69 
70  //Parameters for the total number of SimClusters per layer
71  minTotNsimClsperlay_(pset.getParameter<double>("minTotNsimClsperlay")),
72  maxTotNsimClsperlay_(pset.getParameter<double>("maxTotNsimClsperlay")),
73  nintTotNsimClsperlay_(pset.getParameter<int>("nintTotNsimClsperlay")),
74 
75  //Parameters for the total number of layer clusters per layer
76  minTotNClsperlay_(pset.getParameter<double>("minTotNClsperlay")),
77  maxTotNClsperlay_(pset.getParameter<double>("maxTotNClsperlay")),
78  nintTotNClsperlay_(pset.getParameter<int>("nintTotNClsperlay")),
79 
80  //Parameters for the energy clustered by layer clusters per layer (fraction over CaloParticless)
81  minEneClperlay_(pset.getParameter<double>("minEneClperlay")),
82  maxEneClperlay_(pset.getParameter<double>("maxEneClperlay")),
83  nintEneClperlay_(pset.getParameter<int>("nintEneClperlay")),
84 
85  //Parameters for the score both for:
86  //1. calo particle to layer clusters association per layer
87  //2. layer cluster to calo particles association per layer
88  minScore_(pset.getParameter<double>("minScore")),
89  maxScore_(pset.getParameter<double>("maxScore")),
90  nintScore_(pset.getParameter<int>("nintScore")),
91 
92  //Parameters for shared energy fraction. That is:
93  //1. Fraction of each of the layer clusters energy related to a
94  //calo particle over that calo particle's energy.
95  //2. Fraction of each of the calo particles energy
96  //related to a layer cluster over that layer cluster's energy.
97  minSharedEneFrac_(pset.getParameter<double>("minSharedEneFrac")),
98  maxSharedEneFrac_(pset.getParameter<double>("maxSharedEneFrac")),
99  nintSharedEneFrac_(pset.getParameter<int>("nintSharedEneFrac")),
100  minTSTSharedEneFracEfficiency_(pset.getParameter<double>("minTSTSharedEneFracEfficiency")),
101 
102  //Same as above for Tracksters
103  minTSTSharedEneFrac_(pset.getParameter<double>("minTSTSharedEneFrac")),
104  maxTSTSharedEneFrac_(pset.getParameter<double>("maxTSTSharedEneFrac")),
105  nintTSTSharedEneFrac_(pset.getParameter<int>("nintTSTSharedEneFrac")),
106 
107  //Parameters for the total number of SimClusters per thickness
108  minTotNsimClsperthick_(pset.getParameter<double>("minTotNsimClsperthick")),
109  maxTotNsimClsperthick_(pset.getParameter<double>("maxTotNsimClsperthick")),
110  nintTotNsimClsperthick_(pset.getParameter<int>("nintTotNsimClsperthick")),
111 
112  //Parameters for the total number of layer clusters per thickness
113  minTotNClsperthick_(pset.getParameter<double>("minTotNClsperthick")),
114  maxTotNClsperthick_(pset.getParameter<double>("maxTotNClsperthick")),
115  nintTotNClsperthick_(pset.getParameter<int>("nintTotNClsperthick")),
116 
117  //Parameters for the total number of cells per per thickness per layer
118  minTotNcellsperthickperlayer_(pset.getParameter<double>("minTotNcellsperthickperlayer")),
119  maxTotNcellsperthickperlayer_(pset.getParameter<double>("maxTotNcellsperthickperlayer")),
120  nintTotNcellsperthickperlayer_(pset.getParameter<int>("nintTotNcellsperthickperlayer")),
121 
122  //Parameters for the distance of cluster cells to seed cell per thickness per layer
123  minDisToSeedperthickperlayer_(pset.getParameter<double>("minDisToSeedperthickperlayer")),
124  maxDisToSeedperthickperlayer_(pset.getParameter<double>("maxDisToSeedperthickperlayer")),
125  nintDisToSeedperthickperlayer_(pset.getParameter<int>("nintDisToSeedperthickperlayer")),
126 
127  //Parameters for the energy weighted distance of cluster cells to seed cell per thickness per layer
128  minDisToSeedperthickperlayerenewei_(pset.getParameter<double>("minDisToSeedperthickperlayerenewei")),
129  maxDisToSeedperthickperlayerenewei_(pset.getParameter<double>("maxDisToSeedperthickperlayerenewei")),
130  nintDisToSeedperthickperlayerenewei_(pset.getParameter<int>("nintDisToSeedperthickperlayerenewei")),
131 
132  //Parameters for the distance of cluster cells to max cell per thickness per layer
133  minDisToMaxperthickperlayer_(pset.getParameter<double>("minDisToMaxperthickperlayer")),
134  maxDisToMaxperthickperlayer_(pset.getParameter<double>("maxDisToMaxperthickperlayer")),
135  nintDisToMaxperthickperlayer_(pset.getParameter<int>("nintDisToMaxperthickperlayer")),
136 
137  //Parameters for the energy weighted distance of cluster cells to max cell per thickness per layer
138  minDisToMaxperthickperlayerenewei_(pset.getParameter<double>("minDisToMaxperthickperlayerenewei")),
139  maxDisToMaxperthickperlayerenewei_(pset.getParameter<double>("maxDisToMaxperthickperlayerenewei")),
140  nintDisToMaxperthickperlayerenewei_(pset.getParameter<int>("nintDisToMaxperthickperlayerenewei")),
141 
142  //Parameters for the distance of seed cell to max cell per thickness per layer
143  minDisSeedToMaxperthickperlayer_(pset.getParameter<double>("minDisSeedToMaxperthickperlayer")),
144  maxDisSeedToMaxperthickperlayer_(pset.getParameter<double>("maxDisSeedToMaxperthickperlayer")),
145  nintDisSeedToMaxperthickperlayer_(pset.getParameter<int>("nintDisSeedToMaxperthickperlayer")),
146 
147  //Parameters for the energy of a cluster per thickness per layer
148  minClEneperthickperlayer_(pset.getParameter<double>("minClEneperthickperlayer")),
149  maxClEneperthickperlayer_(pset.getParameter<double>("maxClEneperthickperlayer")),
150  nintClEneperthickperlayer_(pset.getParameter<int>("nintClEneperthickperlayer")),
151 
152  //Parameters for the energy density of cluster cells per thickness
153  minCellsEneDensperthick_(pset.getParameter<double>("minCellsEneDensperthick")),
154  maxCellsEneDensperthick_(pset.getParameter<double>("maxCellsEneDensperthick")),
155  nintCellsEneDensperthick_(pset.getParameter<int>("nintCellsEneDensperthick")),
156 
157  //Parameters for the total number of Tracksters per event
158  // Always treat one event as two events, one in +z one in -z
159  minTotNTSTs_(pset.getParameter<double>("minTotNTSTs")),
160  maxTotNTSTs_(pset.getParameter<double>("maxTotNTSTs")),
161  nintTotNTSTs_(pset.getParameter<int>("nintTotNTSTs")),
162 
163  //Parameters for the total number of layer clusters in Trackster
164  minTotNClsinTSTs_(pset.getParameter<double>("minTotNClsinTSTs")),
165  maxTotNClsinTSTs_(pset.getParameter<double>("maxTotNClsinTSTs")),
166  nintTotNClsinTSTs_(pset.getParameter<int>("nintTotNClsinTSTs")),
167 
168  //Parameters for the total number of layer clusters in Trackster per layer
169  minTotNClsinTSTsperlayer_(pset.getParameter<double>("minTotNClsinTSTsperlayer")),
170  maxTotNClsinTSTsperlayer_(pset.getParameter<double>("maxTotNClsinTSTsperlayer")),
171  nintTotNClsinTSTsperlayer_(pset.getParameter<int>("nintTotNClsinTSTsperlayer")),
172 
173  //Parameters for the multiplicity of layer clusters in Trackster
174  minMplofLCs_(pset.getParameter<double>("minMplofLCs")),
175  maxMplofLCs_(pset.getParameter<double>("maxMplofLCs")),
176  nintMplofLCs_(pset.getParameter<int>("nintMplofLCs")),
177 
178  //Parameters for cluster size
179  minSizeCLsinTSTs_(pset.getParameter<double>("minSizeCLsinTSTs")),
180  maxSizeCLsinTSTs_(pset.getParameter<double>("maxSizeCLsinTSTs")),
181  nintSizeCLsinTSTs_(pset.getParameter<int>("nintSizeCLsinTSTs")),
182 
183  //Parameters for the energy of a cluster per thickness per layer
184  minClEnepermultiplicity_(pset.getParameter<double>("minClEnepermultiplicity")),
185  maxClEnepermultiplicity_(pset.getParameter<double>("maxClEnepermultiplicity")),
186  nintClEnepermultiplicity_(pset.getParameter<int>("nintClEnepermultiplicity")),
187 
188  //parameters for x
189  minX_(pset.getParameter<double>("minX")),
190  maxX_(pset.getParameter<double>("maxX")),
191  nintX_(pset.getParameter<int>("nintX")),
192 
193  //parameters for y
194  minY_(pset.getParameter<double>("minY")),
195  maxY_(pset.getParameter<double>("maxY")),
196  nintY_(pset.getParameter<int>("nintY")),
197 
198  //parameters for z
199  minZ_(pset.getParameter<double>("minZ")),
200  maxZ_(pset.getParameter<double>("maxZ")),
201  nintZ_(pset.getParameter<int>("nintZ")) {}
202 
203 HGVHistoProducerAlgo::~HGVHistoProducerAlgo() {}
204 
205 void HGVHistoProducerAlgo::bookInfo(DQMStore::IBooker& ibook, Histograms& histograms) {
206  histograms.lastLayerEEzm = ibook.bookInt("lastLayerEEzm");
207  histograms.lastLayerFHzm = ibook.bookInt("lastLayerFHzm");
208  histograms.maxlayerzm = ibook.bookInt("maxlayerzm");
209  histograms.lastLayerEEzp = ibook.bookInt("lastLayerEEzp");
210  histograms.lastLayerFHzp = ibook.bookInt("lastLayerFHzp");
211  histograms.maxlayerzp = ibook.bookInt("maxlayerzp");
212 }
213 
214 void HGVHistoProducerAlgo::bookCaloParticleHistos(DQMStore::IBooker& ibook,
215  Histograms& histograms,
216  int pdgid,
217  unsigned int layers) {
218  histograms.h_caloparticle_eta[pdgid] =
219  ibook.book1D("N of caloparticle vs eta", "N of caloParticles vs eta", nintEta_, minEta_, maxEta_);
220  histograms.h_caloparticle_eta_Zorigin[pdgid] =
221  ibook.book2D("Eta vs Zorigin", "Eta vs Zorigin", nintEta_, minEta_, maxEta_, nintZpos_, minZpos_, maxZpos_);
222 
223  histograms.h_caloparticle_energy[pdgid] =
224  ibook.book1D("Energy", "Energy of CaloParticles; Energy [GeV]", nintEne_, minEne_, maxEne_);
225  histograms.h_caloparticle_pt[pdgid] = ibook.book1D("Pt", "Pt of CaloParticles", nintPt_, minPt_, maxPt_);
226  histograms.h_caloparticle_phi[pdgid] = ibook.book1D("Phi", "Phi of CaloParticles", nintPhi_, minPhi_, maxPhi_);
227  histograms.h_caloparticle_selfenergy[pdgid] =
228  ibook.book1D("SelfEnergy", "Total Energy of Hits in Sim Clusters (matched)", nintEne_, minEne_, maxEne_);
229  histograms.h_caloparticle_energyDifference[pdgid] =
230  ibook.book1D("EnergyDifference", "(Energy-SelfEnergy)/Energy", 300., -5., 1.);
231 
232  histograms.h_caloparticle_nSimClusters[pdgid] =
233  ibook.book1D("Num Sim Clusters", "Num Sim Clusters in CaloParticles", 100, 0., 100.);
234  histograms.h_caloparticle_nHitsInSimClusters[pdgid] =
235  ibook.book1D("Num Hits in Sim Clusters", "Num Hits in Sim Clusters in CaloParticles", 1000, 0., 1000.);
236  histograms.h_caloparticle_nHitsInSimClusters_matchedtoRecHit[pdgid] = ibook.book1D(
237  "Num Rec-matched Hits in Sim Clusters", "Num Hits in Sim Clusters (matched) in CaloParticles", 1000, 0., 1000.);
238 
239  histograms.h_caloparticle_nHits_matched_energy[pdgid] =
240  ibook.book1D("Energy of Rec-matched Hits", "Energy of Hits in Sim Clusters (matched)", 100, 0., 10.);
241  histograms.h_caloparticle_nHits_matched_energy_layer[pdgid] =
242  ibook.book2D("Energy of Rec-matched Hits vs layer",
243  "Energy of Hits in Sim Clusters (matched) vs layer",
244  2 * layers,
245  0.,
246  (float)2 * layers,
247  100,
248  0.,
249  10.);
250  histograms.h_caloparticle_nHits_matched_energy_layer_1SimCl[pdgid] =
251  ibook.book2D("Energy of Rec-matched Hits vs layer (1SC)",
252  "Energy of Hits only 1 Sim Clusters (matched) vs layer",
253  2 * layers,
254  0.,
255  (float)2 * layers,
256  100,
257  0.,
258  10.);
259  histograms.h_caloparticle_sum_energy_layer[pdgid] =
260  ibook.book2D("Rec-matched Hits Sum Energy vs layer",
261  "Rescaled Sum Energy of Hits in Sim Clusters (matched) vs layer",
262  2 * layers,
263  0.,
264  (float)2 * layers,
265  110,
266  0.,
267  110.);
268  histograms.h_caloparticle_fractions[pdgid] =
269  ibook.book2D("HitFractions", "Hit fractions;Hit fraction;E_{hit}^{2} fraction", 101, 0, 1.01, 100, 0, 1);
270  histograms.h_caloparticle_fractions_weight[pdgid] = ibook.book2D(
271  "HitFractions_weighted", "Hit fractions weighted;Hit fraction;E_{hit}^{2} fraction", 101, 0, 1.01, 100, 0, 1);
272 
273  histograms.h_caloparticle_firstlayer[pdgid] =
274  ibook.book1D("First Layer", "First layer of the CaloParticles", 2 * layers, 0., (float)2 * layers);
275  histograms.h_caloparticle_lastlayer[pdgid] =
276  ibook.book1D("Last Layer", "Last layer of the CaloParticles", 2 * layers, 0., (float)2 * layers);
277  histograms.h_caloparticle_layersnum[pdgid] =
278  ibook.book1D("Number of Layers", "Number of layers of the CaloParticles", 2 * layers, 0., (float)2 * layers);
279  histograms.h_caloparticle_firstlayer_matchedtoRecHit[pdgid] = ibook.book1D(
280  "First Layer (rec-matched hit)", "First layer of the CaloParticles (matched)", 2 * layers, 0., (float)2 * layers);
281  histograms.h_caloparticle_lastlayer_matchedtoRecHit[pdgid] = ibook.book1D(
282  "Last Layer (rec-matched hit)", "Last layer of the CaloParticles (matched)", 2 * layers, 0., (float)2 * layers);
283  histograms.h_caloparticle_layersnum_matchedtoRecHit[pdgid] =
284  ibook.book1D("Number of Layers (rec-matched hit)",
285  "Number of layers of the CaloParticles (matched)",
286  2 * layers,
287  0.,
288  (float)2 * layers);
289 }
290 
291 void HGVHistoProducerAlgo::bookSimClusterHistos(DQMStore::IBooker& ibook,
292  Histograms& histograms,
293  unsigned int layers,
294  std::vector<int> thicknesses) {
295  //---------------------------------------------------------------------------------------------------------------------------
296  for (unsigned ilayer = 0; ilayer < 2 * layers; ++ilayer) {
297  auto istr1 = std::to_string(ilayer);
298  while (istr1.size() < 2) {
299  istr1.insert(0, "0");
300  }
301  // Make a mapping to the regural layer naming plus z- or z+ for convenience
302  std::string istr2 = "";
303  // first with the -z endcap
304  if (ilayer < layers) {
305  istr2 = std::to_string(ilayer + 1) + " in z-";
306  } else { // then for the +z
307  istr2 = std::to_string(ilayer - (layers - 1)) + " in z+";
308  }
309  histograms.h_simclusternum_perlayer[ilayer] = ibook.book1D("totsimclusternum_layer_" + istr1,
310  "total number of SimClusters for layer " + istr2,
311  nintTotNsimClsperlay_,
312  minTotNsimClsperlay_,
313  maxTotNsimClsperlay_);
314 
315  } //end of loop over layers
316  //---------------------------------------------------------------------------------------------------------------------------
317  for (std::vector<int>::iterator it = thicknesses.begin(); it != thicknesses.end(); ++it) {
318  auto istr = std::to_string(*it);
319  histograms.h_simclusternum_perthick[(*it)] = ibook.book1D("totsimclusternum_thick_" + istr,
320  "total number of simclusters for thickness " + istr,
321  nintTotNsimClsperthick_,
322  minTotNsimClsperthick_,
323  maxTotNsimClsperthick_);
324  } //end of loop over thicknesses
325 
326  //---------------------------------------------------------------------------------------------------------------------------
327  //z-
328  histograms.h_mixedhitssimcluster_zminus =
329  ibook.book1D("mixedhitssimcluster_zminus",
330  "N of simclusters that contain hits of more than one kind in z-",
331  nintMixedHitsSimCluster_,
332  minMixedHitsSimCluster_,
333  maxMixedHitsSimCluster_);
334  //z+
335  histograms.h_mixedhitssimcluster_zplus =
336  ibook.book1D("mixedhitssimcluster_zplus",
337  "N of simclusters that contain hits of more than one kind in z+",
338  nintMixedHitsSimCluster_,
339  minMixedHitsSimCluster_,
340  maxMixedHitsSimCluster_);
341 }
342 
343 void HGVHistoProducerAlgo::bookSimClusterAssociationHistos(DQMStore::IBooker& ibook,
344  Histograms& histograms,
345  unsigned int layers,
346  std::vector<int> thicknesses) {
347  std::unordered_map<int, dqm::reco::MonitorElement*> denom_layercl_in_simcl_eta_perlayer;
348  denom_layercl_in_simcl_eta_perlayer.clear();
349  std::unordered_map<int, dqm::reco::MonitorElement*> denom_layercl_in_simcl_phi_perlayer;
350  denom_layercl_in_simcl_phi_perlayer.clear();
351  std::unordered_map<int, dqm::reco::MonitorElement*> score_layercl2simcluster_perlayer;
352  score_layercl2simcluster_perlayer.clear();
353  std::unordered_map<int, dqm::reco::MonitorElement*> sharedenergy_layercl2simcluster_perlayer;
354  sharedenergy_layercl2simcluster_perlayer.clear();
355  std::unordered_map<int, dqm::reco::MonitorElement*> energy_vs_score_layercl2simcluster_perlayer;
356  energy_vs_score_layercl2simcluster_perlayer.clear();
357  std::unordered_map<int, dqm::reco::MonitorElement*> num_layercl_in_simcl_eta_perlayer;
358  num_layercl_in_simcl_eta_perlayer.clear();
359  std::unordered_map<int, dqm::reco::MonitorElement*> num_layercl_in_simcl_phi_perlayer;
360  num_layercl_in_simcl_phi_perlayer.clear();
361  std::unordered_map<int, dqm::reco::MonitorElement*> numMerge_layercl_in_simcl_eta_perlayer;
362  numMerge_layercl_in_simcl_eta_perlayer.clear();
363  std::unordered_map<int, dqm::reco::MonitorElement*> numMerge_layercl_in_simcl_phi_perlayer;
364  numMerge_layercl_in_simcl_phi_perlayer.clear();
365  std::unordered_map<int, dqm::reco::MonitorElement*> sharedenergy_layercl2simcluster_vs_eta_perlayer;
366  sharedenergy_layercl2simcluster_vs_eta_perlayer.clear();
367  std::unordered_map<int, dqm::reco::MonitorElement*> sharedenergy_layercl2simcluster_vs_phi_perlayer;
368  sharedenergy_layercl2simcluster_vs_phi_perlayer.clear();
369  std::unordered_map<int, dqm::reco::MonitorElement*> denom_simcluster_eta_perlayer;
370  denom_simcluster_eta_perlayer.clear();
371  std::unordered_map<int, dqm::reco::MonitorElement*> denom_simcluster_phi_perlayer;
372  denom_simcluster_phi_perlayer.clear();
373  std::unordered_map<int, dqm::reco::MonitorElement*> score_simcluster2layercl_perlayer;
374  score_simcluster2layercl_perlayer.clear();
375  std::unordered_map<int, dqm::reco::MonitorElement*> sharedenergy_simcluster2layercl_perlayer;
376  sharedenergy_simcluster2layercl_perlayer.clear();
377  std::unordered_map<int, dqm::reco::MonitorElement*> energy_vs_score_simcluster2layercl_perlayer;
378  energy_vs_score_simcluster2layercl_perlayer.clear();
379  std::unordered_map<int, dqm::reco::MonitorElement*> num_simcluster_eta_perlayer;
380  num_simcluster_eta_perlayer.clear();
381  std::unordered_map<int, dqm::reco::MonitorElement*> num_simcluster_phi_perlayer;
382  num_simcluster_phi_perlayer.clear();
383  std::unordered_map<int, dqm::reco::MonitorElement*> numDup_simcluster_eta_perlayer;
384  numDup_simcluster_eta_perlayer.clear();
385  std::unordered_map<int, dqm::reco::MonitorElement*> numDup_simcluster_phi_perlayer;
386  numDup_simcluster_phi_perlayer.clear();
387  std::unordered_map<int, dqm::reco::MonitorElement*> sharedenergy_simcluster2layercl_vs_eta_perlayer;
388  sharedenergy_simcluster2layercl_vs_eta_perlayer.clear();
389  std::unordered_map<int, dqm::reco::MonitorElement*> sharedenergy_simcluster2layercl_vs_phi_perlayer;
390  sharedenergy_simcluster2layercl_vs_phi_perlayer.clear();
391 
392  //---------------------------------------------------------------------------------------------------------------------------
393  for (unsigned ilayer = 0; ilayer < 2 * layers; ++ilayer) {
394  auto istr1 = std::to_string(ilayer);
395  while (istr1.size() < 2) {
396  istr1.insert(0, "0");
397  }
398  // Make a mapping to the regural layer naming plus z- or z+ for convenience
399  std::string istr2 = "";
400  // first with the -z endcap
401  if (ilayer < layers) {
402  istr2 = std::to_string(ilayer + 1) + " in z-";
403  } else { // then for the +z
404  istr2 = std::to_string(ilayer - (layers - 1)) + " in z+";
405  }
406  //-------------------------------------------------------------------------------------------------------------------------
407  denom_layercl_in_simcl_eta_perlayer[ilayer] =
408  ibook.book1D("Denom_LayerCluster_in_SimCluster_Eta_perlayer" + istr1,
409  "Denom LayerCluster in SimCluster Eta per Layer Cluster for layer " + istr2,
410  nintEta_,
411  minEta_,
412  maxEta_);
413  //-------------------------------------------------------------------------------------------------------------------------
414  denom_layercl_in_simcl_phi_perlayer[ilayer] =
415  ibook.book1D("Denom_LayerCluster_in_SimCluster_Phi_perlayer" + istr1,
416  "Denom LayerCluster in SimCluster Phi per Layer Cluster for layer " + istr2,
417  nintPhi_,
418  minPhi_,
419  maxPhi_);
420  //-------------------------------------------------------------------------------------------------------------------------
421  score_layercl2simcluster_perlayer[ilayer] = ibook.book1D("Score_layercl2simcluster_perlayer" + istr1,
422  "Score of Layer Cluster per SimCluster for layer " + istr2,
423  nintScore_,
424  minScore_,
425  maxScore_);
426  //-------------------------------------------------------------------------------------------------------------------------
427  score_simcluster2layercl_perlayer[ilayer] = ibook.book1D("Score_simcluster2layercl_perlayer" + istr1,
428  "Score of SimCluster per Layer Cluster for layer " + istr2,
429  nintScore_,
430  minScore_,
431  maxScore_);
432  //-------------------------------------------------------------------------------------------------------------------------
433  energy_vs_score_simcluster2layercl_perlayer[ilayer] =
434  ibook.book2D("Energy_vs_Score_simcluster2layer_perlayer" + istr1,
435  "Energy vs Score of SimCluster per Layer Cluster for layer " + istr2,
436  nintScore_,
437  minScore_,
438  maxScore_,
439  nintSharedEneFrac_,
440  minSharedEneFrac_,
441  maxSharedEneFrac_);
442  //-------------------------------------------------------------------------------------------------------------------------
443  energy_vs_score_layercl2simcluster_perlayer[ilayer] =
444  ibook.book2D("Energy_vs_Score_layer2simcluster_perlayer" + istr1,
445  "Energy vs Score of Layer Cluster per SimCluster for layer " + istr2,
446  nintScore_,
447  minScore_,
448  maxScore_,
449  nintSharedEneFrac_,
450  minSharedEneFrac_,
451  maxSharedEneFrac_);
452  //-------------------------------------------------------------------------------------------------------------------------
453  sharedenergy_simcluster2layercl_perlayer[ilayer] =
454  ibook.book1D("SharedEnergy_simcluster2layercl_perlayer" + istr1,
455  "Shared Energy of SimCluster per Layer Cluster for layer " + istr2,
456  nintSharedEneFrac_,
457  minSharedEneFrac_,
458  maxSharedEneFrac_);
459  //-------------------------------------------------------------------------------------------------------------------------
460  sharedenergy_simcluster2layercl_vs_eta_perlayer[ilayer] =
461  ibook.bookProfile("SharedEnergy_simcluster2layercl_vs_eta_perlayer" + istr1,
462  "Shared Energy of SimCluster vs #eta per best Layer Cluster for layer " + istr2,
463  nintEta_,
464  minEta_,
465  maxEta_,
466  minSharedEneFrac_,
467  maxSharedEneFrac_);
468  //-------------------------------------------------------------------------------------------------------------------------
469  sharedenergy_simcluster2layercl_vs_phi_perlayer[ilayer] =
470  ibook.bookProfile("SharedEnergy_simcluster2layercl_vs_phi_perlayer" + istr1,
471  "Shared Energy of SimCluster vs #phi per best Layer Cluster for layer " + istr2,
472  nintPhi_,
473  minPhi_,
474  maxPhi_,
475  minSharedEneFrac_,
476  maxSharedEneFrac_);
477  //-------------------------------------------------------------------------------------------------------------------------
478  sharedenergy_layercl2simcluster_perlayer[ilayer] =
479  ibook.book1D("SharedEnergy_layercluster2simcluster_perlayer" + istr1,
480  "Shared Energy of Layer Cluster per SimCluster for layer " + istr2,
481  nintSharedEneFrac_,
482  minSharedEneFrac_,
483  maxSharedEneFrac_);
484  //-------------------------------------------------------------------------------------------------------------------------
485  sharedenergy_layercl2simcluster_vs_eta_perlayer[ilayer] =
486  ibook.bookProfile("SharedEnergy_layercl2simcluster_vs_eta_perlayer" + istr1,
487  "Shared Energy of LayerCluster vs #eta per best SimCluster for layer " + istr2,
488  nintEta_,
489  minEta_,
490  maxEta_,
491  minSharedEneFrac_,
492  maxSharedEneFrac_);
493  //-------------------------------------------------------------------------------------------------------------------------
494  sharedenergy_layercl2simcluster_vs_phi_perlayer[ilayer] =
495  ibook.bookProfile("SharedEnergy_layercl2simcluster_vs_phi_perlayer" + istr1,
496  "Shared Energy of LayerCluster vs #phi per best SimCluster for layer " + istr2,
497  nintPhi_,
498  minPhi_,
499  maxPhi_,
500  minSharedEneFrac_,
501  maxSharedEneFrac_);
502  //-------------------------------------------------------------------------------------------------------------------------
503  num_simcluster_eta_perlayer[ilayer] = ibook.book1D("Num_SimCluster_Eta_perlayer" + istr1,
504  "Num SimCluster Eta per Layer Cluster for layer " + istr2,
505  nintEta_,
506  minEta_,
507  maxEta_);
508  //-------------------------------------------------------------------------------------------------------------------------
509  numDup_simcluster_eta_perlayer[ilayer] =
510  ibook.book1D("NumDup_SimCluster_Eta_perlayer" + istr1,
511  "Num Duplicate SimCluster Eta per Layer Cluster for layer " + istr2,
512  nintEta_,
513  minEta_,
514  maxEta_);
515  //-------------------------------------------------------------------------------------------------------------------------
516  denom_simcluster_eta_perlayer[ilayer] = ibook.book1D("Denom_SimCluster_Eta_perlayer" + istr1,
517  "Denom SimCluster Eta per Layer Cluster for layer " + istr2,
518  nintEta_,
519  minEta_,
520  maxEta_);
521  //-------------------------------------------------------------------------------------------------------------------------
522  num_simcluster_phi_perlayer[ilayer] = ibook.book1D("Num_SimCluster_Phi_perlayer" + istr1,
523  "Num SimCluster Phi per Layer Cluster for layer " + istr2,
524  nintPhi_,
525  minPhi_,
526  maxPhi_);
527  //-------------------------------------------------------------------------------------------------------------------------
528  numDup_simcluster_phi_perlayer[ilayer] =
529  ibook.book1D("NumDup_SimCluster_Phi_perlayer" + istr1,
530  "Num Duplicate SimCluster Phi per Layer Cluster for layer " + istr2,
531  nintPhi_,
532  minPhi_,
533  maxPhi_);
534  //-------------------------------------------------------------------------------------------------------------------------
535  denom_simcluster_phi_perlayer[ilayer] = ibook.book1D("Denom_SimCluster_Phi_perlayer" + istr1,
536  "Denom SimCluster Phi per Layer Cluster for layer " + istr2,
537  nintPhi_,
538  minPhi_,
539  maxPhi_);
540  //-------------------------------------------------------------------------------------------------------------------------
541  num_layercl_in_simcl_eta_perlayer[ilayer] =
542  ibook.book1D("Num_LayerCluster_in_SimCluster_Eta_perlayer" + istr1,
543  "Num LayerCluster Eta per Layer Cluster in SimCluster for layer " + istr2,
544  nintEta_,
545  minEta_,
546  maxEta_);
547  //-------------------------------------------------------------------------------------------------------------------------
548  numMerge_layercl_in_simcl_eta_perlayer[ilayer] =
549  ibook.book1D("NumMerge_LayerCluster_in_SimCluster_Eta_perlayer" + istr1,
550  "Num Merge LayerCluster Eta per Layer Cluster in SimCluster for layer " + istr2,
551  nintEta_,
552  minEta_,
553  maxEta_);
554  //-------------------------------------------------------------------------------------------------------------------------
555  num_layercl_in_simcl_phi_perlayer[ilayer] =
556  ibook.book1D("Num_LayerCluster_in_SimCluster_Phi_perlayer" + istr1,
557  "Num LayerCluster Phi per Layer Cluster in SimCluster for layer " + istr2,
558  nintPhi_,
559  minPhi_,
560  maxPhi_);
561  //-------------------------------------------------------------------------------------------------------------------------
562  numMerge_layercl_in_simcl_phi_perlayer[ilayer] =
563  ibook.book1D("NumMerge_LayerCluster_in_SimCluster_Phi_perlayer" + istr1,
564  "Num Merge LayerCluster Phi per Layer Cluster in SimCluster for layer " + istr2,
565  nintPhi_,
566  minPhi_,
567  maxPhi_);
568 
569  } //end of loop over layers
570 
571  histograms.h_denom_layercl_in_simcl_eta_perlayer.push_back(std::move(denom_layercl_in_simcl_eta_perlayer));
572  histograms.h_denom_layercl_in_simcl_phi_perlayer.push_back(std::move(denom_layercl_in_simcl_phi_perlayer));
573  histograms.h_score_layercl2simcluster_perlayer.push_back(std::move(score_layercl2simcluster_perlayer));
574  histograms.h_sharedenergy_layercl2simcluster_perlayer.push_back(std::move(sharedenergy_layercl2simcluster_perlayer));
575  histograms.h_energy_vs_score_layercl2simcluster_perlayer.push_back(
576  std::move(energy_vs_score_layercl2simcluster_perlayer));
577  histograms.h_num_layercl_in_simcl_eta_perlayer.push_back(std::move(num_layercl_in_simcl_eta_perlayer));
578  histograms.h_num_layercl_in_simcl_phi_perlayer.push_back(std::move(num_layercl_in_simcl_phi_perlayer));
579  histograms.h_numMerge_layercl_in_simcl_eta_perlayer.push_back(std::move(numMerge_layercl_in_simcl_eta_perlayer));
580  histograms.h_numMerge_layercl_in_simcl_phi_perlayer.push_back(std::move(numMerge_layercl_in_simcl_phi_perlayer));
581  histograms.h_sharedenergy_layercl2simcluster_vs_eta_perlayer.push_back(
582  std::move(sharedenergy_layercl2simcluster_vs_eta_perlayer));
583  histograms.h_sharedenergy_layercl2simcluster_vs_phi_perlayer.push_back(
584  std::move(sharedenergy_layercl2simcluster_vs_phi_perlayer));
585  histograms.h_denom_simcluster_eta_perlayer.push_back(std::move(denom_simcluster_eta_perlayer));
586  histograms.h_denom_simcluster_phi_perlayer.push_back(std::move(denom_simcluster_phi_perlayer));
587  histograms.h_score_simcluster2layercl_perlayer.push_back(std::move(score_simcluster2layercl_perlayer));
588  histograms.h_sharedenergy_simcluster2layercl_perlayer.push_back(std::move(sharedenergy_simcluster2layercl_perlayer));
589  histograms.h_energy_vs_score_simcluster2layercl_perlayer.push_back(
590  std::move(energy_vs_score_simcluster2layercl_perlayer));
591  histograms.h_num_simcluster_eta_perlayer.push_back(std::move(num_simcluster_eta_perlayer));
592  histograms.h_num_simcluster_phi_perlayer.push_back(std::move(num_simcluster_phi_perlayer));
593  histograms.h_numDup_simcluster_eta_perlayer.push_back(std::move(numDup_simcluster_eta_perlayer));
594  histograms.h_numDup_simcluster_phi_perlayer.push_back(std::move(numDup_simcluster_phi_perlayer));
595  histograms.h_sharedenergy_simcluster2layercl_vs_eta_perlayer.push_back(
596  std::move(sharedenergy_simcluster2layercl_vs_eta_perlayer));
597  histograms.h_sharedenergy_simcluster2layercl_vs_phi_perlayer.push_back(
598  std::move(sharedenergy_simcluster2layercl_vs_phi_perlayer));
599 }
600 void HGVHistoProducerAlgo::bookClusterHistos_ClusterLevel(DQMStore::IBooker& ibook,
601  Histograms& histograms,
602  unsigned int layers,
603  std::vector<int> thicknesses,
604  std::string pathtomatbudfile) {
605  //---------------------------------------------------------------------------------------------------------------------------
606  histograms.h_cluster_eta.push_back(
607  ibook.book1D("num_reco_cluster_eta", "N of reco clusters vs eta", nintEta_, minEta_, maxEta_));
608 
609  //---------------------------------------------------------------------------------------------------------------------------
610  //z-
611  histograms.h_mixedhitscluster_zminus.push_back(
612  ibook.book1D("mixedhitscluster_zminus",
613  "N of reco clusters that contain hits of more than one kind in z-",
614  nintMixedHitsCluster_,
615  minMixedHitsCluster_,
616  maxMixedHitsCluster_));
617  //z+
618  histograms.h_mixedhitscluster_zplus.push_back(
619  ibook.book1D("mixedhitscluster_zplus",
620  "N of reco clusters that contain hits of more than one kind in z+",
621  nintMixedHitsCluster_,
622  minMixedHitsCluster_,
623  maxMixedHitsCluster_));
624 
625  //---------------------------------------------------------------------------------------------------------------------------
626  //z-
627  histograms.h_energyclustered_zminus.push_back(
628  ibook.book1D("energyclustered_zminus",
629  "percent of total energy clustered by all layer clusters over CaloParticless energy in z-",
630  nintEneCl_,
631  minEneCl_,
632  maxEneCl_));
633  //z+
634  histograms.h_energyclustered_zplus.push_back(
635  ibook.book1D("energyclustered_zplus",
636  "percent of total energy clustered by all layer clusters over CaloParticless energy in z+",
637  nintEneCl_,
638  minEneCl_,
639  maxEneCl_));
640 
641  //---------------------------------------------------------------------------------------------------------------------------
642  //z-
643  std::string subpathtomat = pathtomatbudfile.substr(pathtomatbudfile.find("Validation"));
644  histograms.h_longdepthbarycentre_zminus.push_back(
645  ibook.book1D("longdepthbarycentre_zminus",
646  "The longitudinal depth barycentre in z- for " + subpathtomat,
647  nintLongDepBary_,
648  minLongDepBary_,
649  maxLongDepBary_));
650  //z+
651  histograms.h_longdepthbarycentre_zplus.push_back(
652  ibook.book1D("longdepthbarycentre_zplus",
653  "The longitudinal depth barycentre in z+ for " + subpathtomat,
654  nintLongDepBary_,
655  minLongDepBary_,
656  maxLongDepBary_));
657 
658  //---------------------------------------------------------------------------------------------------------------------------
659  for (unsigned ilayer = 0; ilayer < 2 * layers; ++ilayer) {
660  auto istr1 = std::to_string(ilayer);
661  while (istr1.size() < 2) {
662  istr1.insert(0, "0");
663  }
664  // Make a mapping to the regural layer naming plus z- or z+ for convenience
665  std::string istr2 = "";
666  // first with the -z endcap
667  if (ilayer < layers) {
668  istr2 = std::to_string(ilayer + 1) + " in z-";
669  } else { // then for the +z
670  istr2 = std::to_string(ilayer - (layers - 1)) + " in z+";
671  }
672  histograms.h_clusternum_perlayer[ilayer] = ibook.book1D("totclusternum_layer_" + istr1,
673  "total number of layer clusters for layer " + istr2,
674  nintTotNClsperlay_,
675  minTotNClsperlay_,
676  maxTotNClsperlay_);
677  histograms.h_energyclustered_perlayer[ilayer] = ibook.book1D(
678  "energyclustered_perlayer" + istr1,
679  "percent of total energy clustered by layer clusters over CaloParticless energy for layer " + istr2,
680  nintEneClperlay_,
681  minEneClperlay_,
682  maxEneClperlay_);
683  }
684 
685  //---------------------------------------------------------------------------------------------------------------------------
686  for (std::vector<int>::iterator it = thicknesses.begin(); it != thicknesses.end(); ++it) {
687  auto istr = std::to_string(*it);
688  histograms.h_clusternum_perthick[(*it)] = ibook.book1D("totclusternum_thick_" + istr,
689  "total number of layer clusters for thickness " + istr,
690  nintTotNClsperthick_,
691  minTotNClsperthick_,
692  maxTotNClsperthick_);
693  }
694  //---------------------------------------------------------------------------------------------------------------------------
695 }
696 
697 void HGVHistoProducerAlgo::bookClusterHistos_LCtoCP_association(DQMStore::IBooker& ibook,
698  Histograms& histograms,
699  unsigned int layers) {
700  //----------------------------------------------------------------------------------------------------------------------------
701  for (unsigned ilayer = 0; ilayer < 2 * layers; ++ilayer) {
702  auto istr1 = std::to_string(ilayer);
703  while (istr1.size() < 2) {
704  istr1.insert(0, "0");
705  }
706  // Make a mapping to the regural layer naming plus z- or z+ for convenience
707  std::string istr2 = "";
708  // first with the -z endcap
709  if (ilayer < layers) {
710  istr2 = std::to_string(ilayer + 1) + " in z-";
711  } else { // then for the +z
712  istr2 = std::to_string(ilayer - (layers - 1)) + " in z+";
713  }
714  histograms.h_score_layercl2caloparticle_perlayer[ilayer] =
715  ibook.book1D("Score_layercl2caloparticle_perlayer" + istr1,
716  "Score of Layer Cluster per CaloParticle for layer " + istr2,
717  nintScore_,
718  minScore_,
719  maxScore_);
720  histograms.h_score_caloparticle2layercl_perlayer[ilayer] =
721  ibook.book1D("Score_caloparticle2layercl_perlayer" + istr1,
722  "Score of CaloParticle per Layer Cluster for layer " + istr2,
723  nintScore_,
724  minScore_,
725  maxScore_);
726  histograms.h_energy_vs_score_caloparticle2layercl_perlayer[ilayer] =
727  ibook.book2D("Energy_vs_Score_caloparticle2layer_perlayer" + istr1,
728  "Energy vs Score of CaloParticle per Layer Cluster for layer " + istr2,
729  nintScore_,
730  minScore_,
731  maxScore_,
732  nintSharedEneFrac_,
733  minSharedEneFrac_,
734  maxSharedEneFrac_);
735  histograms.h_energy_vs_score_layercl2caloparticle_perlayer[ilayer] =
736  ibook.book2D("Energy_vs_Score_layer2caloparticle_perlayer" + istr1,
737  "Energy vs Score of Layer Cluster per CaloParticle Layer for layer " + istr2,
738  nintScore_,
739  minScore_,
740  maxScore_,
741  nintSharedEneFrac_,
742  minSharedEneFrac_,
743  maxSharedEneFrac_);
744  histograms.h_sharedenergy_caloparticle2layercl_perlayer[ilayer] =
745  ibook.book1D("SharedEnergy_caloparticle2layercl_perlayer" + istr1,
746  "Shared Energy of CaloParticle per Layer Cluster for layer " + istr2,
747  nintSharedEneFrac_,
748  minSharedEneFrac_,
749  maxSharedEneFrac_);
750  histograms.h_sharedenergy_caloparticle2layercl_vs_eta_perlayer[ilayer] =
751  ibook.bookProfile("SharedEnergy_caloparticle2layercl_vs_eta_perlayer" + istr1,
752  "Shared Energy of CaloParticle vs #eta per best Layer Cluster for layer " + istr2,
753  nintEta_,
754  minEta_,
755  maxEta_,
756  minSharedEneFrac_,
757  maxSharedEneFrac_);
758  histograms.h_sharedenergy_caloparticle2layercl_vs_phi_perlayer[ilayer] =
759  ibook.bookProfile("SharedEnergy_caloparticle2layercl_vs_phi_perlayer" + istr1,
760  "Shared Energy of CaloParticle vs #phi per best Layer Cluster for layer " + istr2,
761  nintPhi_,
762  minPhi_,
763  maxPhi_,
764  minSharedEneFrac_,
765  maxSharedEneFrac_);
766  histograms.h_sharedenergy_layercl2caloparticle_perlayer[ilayer] =
767  ibook.book1D("SharedEnergy_layercluster2caloparticle_perlayer" + istr1,
768  "Shared Energy of Layer Cluster per Layer Calo Particle for layer " + istr2,
769  nintSharedEneFrac_,
770  minSharedEneFrac_,
771  maxSharedEneFrac_);
772  histograms.h_sharedenergy_layercl2caloparticle_vs_eta_perlayer[ilayer] =
773  ibook.bookProfile("SharedEnergy_layercl2caloparticle_vs_eta_perlayer" + istr1,
774  "Shared Energy of LayerCluster vs #eta per best Calo Particle for layer " + istr2,
775  nintEta_,
776  minEta_,
777  maxEta_,
778  minSharedEneFrac_,
779  maxSharedEneFrac_);
780  histograms.h_sharedenergy_layercl2caloparticle_vs_phi_perlayer[ilayer] =
781  ibook.bookProfile("SharedEnergy_layercl2caloparticle_vs_phi_perlayer" + istr1,
782  "Shared Energy of LayerCluster vs #phi per best Calo Particle for layer " + istr2,
783  nintPhi_,
784  minPhi_,
785  maxPhi_,
786  minSharedEneFrac_,
787  maxSharedEneFrac_);
788  histograms.h_num_caloparticle_eta_perlayer[ilayer] =
789  ibook.book1D("Num_CaloParticle_Eta_perlayer" + istr1,
790  "Num CaloParticle Eta per Layer Cluster for layer " + istr2,
791  nintEta_,
792  minEta_,
793  maxEta_);
794  histograms.h_numDup_caloparticle_eta_perlayer[ilayer] =
795  ibook.book1D("NumDup_CaloParticle_Eta_perlayer" + istr1,
796  "Num Duplicate CaloParticle Eta per Layer Cluster for layer " + istr2,
797  nintEta_,
798  minEta_,
799  maxEta_);
800  histograms.h_denom_caloparticle_eta_perlayer[ilayer] =
801  ibook.book1D("Denom_CaloParticle_Eta_perlayer" + istr1,
802  "Denom CaloParticle Eta per Layer Cluster for layer " + istr2,
803  nintEta_,
804  minEta_,
805  maxEta_);
806  histograms.h_num_caloparticle_phi_perlayer[ilayer] =
807  ibook.book1D("Num_CaloParticle_Phi_perlayer" + istr1,
808  "Num CaloParticle Phi per Layer Cluster for layer " + istr2,
809  nintPhi_,
810  minPhi_,
811  maxPhi_);
812  histograms.h_numDup_caloparticle_phi_perlayer[ilayer] =
813  ibook.book1D("NumDup_CaloParticle_Phi_perlayer" + istr1,
814  "Num Duplicate CaloParticle Phi per Layer Cluster for layer " + istr2,
815  nintPhi_,
816  minPhi_,
817  maxPhi_);
818  histograms.h_denom_caloparticle_phi_perlayer[ilayer] =
819  ibook.book1D("Denom_CaloParticle_Phi_perlayer" + istr1,
820  "Denom CaloParticle Phi per Layer Cluster for layer " + istr2,
821  nintPhi_,
822  minPhi_,
823  maxPhi_);
824  histograms.h_num_layercl_eta_perlayer[ilayer] =
825  ibook.book1D("Num_LayerCluster_Eta_perlayer" + istr1,
826  "Num LayerCluster Eta per Layer Cluster for layer " + istr2,
827  nintEta_,
828  minEta_,
829  maxEta_);
830  histograms.h_numMerge_layercl_eta_perlayer[ilayer] =
831  ibook.book1D("NumMerge_LayerCluster_Eta_perlayer" + istr1,
832  "Num Merge LayerCluster Eta per Layer Cluster for layer " + istr2,
833  nintEta_,
834  minEta_,
835  maxEta_);
836  histograms.h_denom_layercl_eta_perlayer[ilayer] =
837  ibook.book1D("Denom_LayerCluster_Eta_perlayer" + istr1,
838  "Denom LayerCluster Eta per Layer Cluster for layer " + istr2,
839  nintEta_,
840  minEta_,
841  maxEta_);
842  histograms.h_num_layercl_phi_perlayer[ilayer] =
843  ibook.book1D("Num_LayerCluster_Phi_perlayer" + istr1,
844  "Num LayerCluster Phi per Layer Cluster for layer " + istr2,
845  nintPhi_,
846  minPhi_,
847  maxPhi_);
848  histograms.h_numMerge_layercl_phi_perlayer[ilayer] =
849  ibook.book1D("NumMerge_LayerCluster_Phi_perlayer" + istr1,
850  "Num Merge LayerCluster Phi per Layer Cluster for layer " + istr2,
851  nintPhi_,
852  minPhi_,
853  maxPhi_);
854  histograms.h_denom_layercl_phi_perlayer[ilayer] =
855  ibook.book1D("Denom_LayerCluster_Phi_perlayer" + istr1,
856  "Denom LayerCluster Phi per Layer Cluster for layer " + istr2,
857  nintPhi_,
858  minPhi_,
859  maxPhi_);
860  }
861  //---------------------------------------------------------------------------------------------------------------------------
862 }
863 
864 void HGVHistoProducerAlgo::bookClusterHistos_CellLevel(DQMStore::IBooker& ibook,
865  Histograms& histograms,
866  unsigned int layers,
867  std::vector<int> thicknesses) {
868  //----------------------------------------------------------------------------------------------------------------------------
869  for (unsigned ilayer = 0; ilayer < 2 * layers; ++ilayer) {
870  auto istr1 = std::to_string(ilayer);
871  while (istr1.size() < 2) {
872  istr1.insert(0, "0");
873  }
874  // Make a mapping to the regural layer naming plus z- or z+ for convenience
875  std::string istr2 = "";
876  // first with the -z endcap
877  if (ilayer < layers) {
878  istr2 = std::to_string(ilayer + 1) + " in z-";
879  } else { // then for the +z
880  istr2 = std::to_string(ilayer - (layers - 1)) + " in z+";
881  }
882  histograms.h_cellAssociation_perlayer[ilayer] =
883  ibook.book1D("cellAssociation_perlayer" + istr1, "Cell Association for layer " + istr2, 5, -4., 1.);
884  histograms.h_cellAssociation_perlayer[ilayer]->setBinLabel(2, "TN(purity)");
885  histograms.h_cellAssociation_perlayer[ilayer]->setBinLabel(3, "FN(ineff.)");
886  histograms.h_cellAssociation_perlayer[ilayer]->setBinLabel(4, "FP(fake)");
887  histograms.h_cellAssociation_perlayer[ilayer]->setBinLabel(5, "TP(eff.)");
888  }
889  //----------------------------------------------------------------------------------------------------------------------------
890  for (std::vector<int>::iterator it = thicknesses.begin(); it != thicknesses.end(); ++it) {
891  auto istr = std::to_string(*it);
892  histograms.h_cellsenedens_perthick[(*it)] = ibook.book1D("cellsenedens_thick_" + istr,
893  "energy density of cluster cells for thickness " + istr,
894  nintCellsEneDensperthick_,
895  minCellsEneDensperthick_,
896  maxCellsEneDensperthick_);
897  }
898  //----------------------------------------------------------------------------------------------------------------------------
899  //Not all combination exists but should keep them all for cross checking reason.
900  for (std::vector<int>::iterator it = thicknesses.begin(); it != thicknesses.end(); ++it) {
901  for (unsigned ilayer = 0; ilayer < 2 * layers; ++ilayer) {
902  auto istr1 = std::to_string(*it);
903  auto istr2 = std::to_string(ilayer);
904  while (istr2.size() < 2)
905  istr2.insert(0, "0");
906  auto istr = istr1 + "_" + istr2;
907  // Make a mapping to the regural layer naming plus z- or z+ for convenience
908  std::string istr3 = "";
909  // first with the -z endcap
910  if (ilayer < layers) {
911  istr3 = std::to_string(ilayer + 1) + " in z- ";
912  } else { // then for the +z
913  istr3 = std::to_string(ilayer - (layers - 1)) + " in z+ ";
914  }
915  //---
916  histograms.h_cellsnum_perthickperlayer[istr] =
917  ibook.book1D("cellsnum_perthick_perlayer_" + istr,
918  "total number of cells for layer " + istr3 + " for thickness " + istr1,
919  nintTotNcellsperthickperlayer_,
920  minTotNcellsperthickperlayer_,
921  maxTotNcellsperthickperlayer_);
922  //---
923  histograms.h_distancetoseedcell_perthickperlayer[istr] =
924  ibook.book1D("distancetoseedcell_perthickperlayer_" + istr,
925  "distance of cluster cells to seed cell for layer " + istr3 + " for thickness " + istr1,
926  nintDisToSeedperthickperlayer_,
927  minDisToSeedperthickperlayer_,
928  maxDisToSeedperthickperlayer_);
929  //---
930  histograms.h_distancetoseedcell_perthickperlayer_eneweighted[istr] = ibook.book1D(
931  "distancetoseedcell_perthickperlayer_eneweighted_" + istr,
932  "energy weighted distance of cluster cells to seed cell for layer " + istr3 + " for thickness " + istr1,
933  nintDisToSeedperthickperlayerenewei_,
934  minDisToSeedperthickperlayerenewei_,
935  maxDisToSeedperthickperlayerenewei_);
936  //---
937  histograms.h_distancetomaxcell_perthickperlayer[istr] =
938  ibook.book1D("distancetomaxcell_perthickperlayer_" + istr,
939  "distance of cluster cells to max cell for layer " + istr3 + " for thickness " + istr1,
940  nintDisToMaxperthickperlayer_,
941  minDisToMaxperthickperlayer_,
942  maxDisToMaxperthickperlayer_);
943  //---
944  histograms.h_distancetomaxcell_perthickperlayer_eneweighted[istr] = ibook.book1D(
945  "distancetomaxcell_perthickperlayer_eneweighted_" + istr,
946  "energy weighted distance of cluster cells to max cell for layer " + istr3 + " for thickness " + istr1,
947  nintDisToMaxperthickperlayerenewei_,
948  minDisToMaxperthickperlayerenewei_,
949  maxDisToMaxperthickperlayerenewei_);
950  //---
951  histograms.h_distancebetseedandmaxcell_perthickperlayer[istr] =
952  ibook.book1D("distancebetseedandmaxcell_perthickperlayer_" + istr,
953  "distance of seed cell to max cell for layer " + istr3 + " for thickness " + istr1,
954  nintDisSeedToMaxperthickperlayer_,
955  minDisSeedToMaxperthickperlayer_,
956  maxDisSeedToMaxperthickperlayer_);
957  //---
958  histograms.h_distancebetseedandmaxcellvsclusterenergy_perthickperlayer[istr] = ibook.book2D(
959  "distancebetseedandmaxcellvsclusterenergy_perthickperlayer_" + istr,
960  "distance of seed cell to max cell vs cluster energy for layer " + istr3 + " for thickness " + istr1,
961  nintDisSeedToMaxperthickperlayer_,
962  minDisSeedToMaxperthickperlayer_,
963  maxDisSeedToMaxperthickperlayer_,
964  nintClEneperthickperlayer_,
965  minClEneperthickperlayer_,
966  maxClEneperthickperlayer_);
967  }
968  }
969 }
970 //----------------------------------------------------------------------------------------------------------------------------
971 
972 void HGVHistoProducerAlgo::bookTracksterHistos(DQMStore::IBooker& ibook, Histograms& histograms, unsigned int layers) {
973  std::unordered_map<int, dqm::reco::MonitorElement*> clusternum_in_trackster_perlayer;
974  clusternum_in_trackster_perlayer.clear();
975 
976  for (unsigned ilayer = 0; ilayer < 2 * layers; ++ilayer) {
977  auto istr1 = std::to_string(ilayer);
978  while (istr1.size() < 2) {
979  istr1.insert(0, "0");
980  }
981  // Make a mapping to the regural layer naming plus z- or z+ for convenience
982  std::string istr2 = "";
983  // first with the -z endcap
984  if (ilayer < layers) {
985  istr2 = std::to_string(ilayer + 1) + " in z-";
986  } else { // then for the +z
987  istr2 = std::to_string(ilayer - (layers - 1)) + " in z+";
988  }
989 
990  clusternum_in_trackster_perlayer[ilayer] = ibook.book1D("clusternum_in_trackster_perlayer" + istr1,
991  "Number of layer clusters in Trackster for layer " + istr2,
992  nintTotNClsinTSTsperlayer_,
993  minTotNClsinTSTsperlayer_,
994  maxTotNClsinTSTsperlayer_);
995  }
996 
997  histograms.h_clusternum_in_trackster_perlayer.push_back(std::move(clusternum_in_trackster_perlayer));
998 
999  histograms.h_tracksternum.push_back(ibook.book1D(
1000  "tottracksternum", "total number of Tracksters;# of Tracksters", nintTotNTSTs_, minTotNTSTs_, maxTotNTSTs_));
1001 
1002  histograms.h_conttracksternum.push_back(ibook.book1D(
1003  "conttracksternum", "number of Tracksters with 3 contiguous layers", nintTotNTSTs_, minTotNTSTs_, maxTotNTSTs_));
1004 
1005  histograms.h_nonconttracksternum.push_back(ibook.book1D("nonconttracksternum",
1006  "number of Tracksters without 3 contiguous layers",
1007  nintTotNTSTs_,
1008  minTotNTSTs_,
1009  maxTotNTSTs_));
1010 
1011  histograms.h_clusternum_in_trackster.push_back(
1012  ibook.book1D("clusternum_in_trackster",
1013  "total number of layer clusters in Trackster;# of LayerClusters",
1014  nintTotNClsinTSTs_,
1015  minTotNClsinTSTs_,
1016  maxTotNClsinTSTs_));
1017 
1018  histograms.h_clusternum_in_trackster_vs_layer.push_back(ibook.bookProfile(
1019  "clusternum_in_trackster_vs_layer",
1020  "Profile of 2d layer clusters in Trackster vs layer number;layer number;<2D LayerClusters in Trackster>",
1021  2 * layers,
1022  0.,
1023  2. * layers,
1024  minTotNClsinTSTsperlayer_,
1025  maxTotNClsinTSTsperlayer_));
1026 
1027  histograms.h_multiplicityOfLCinTST.push_back(
1028  ibook.book2D("multiplicityOfLCinTST",
1029  "Multiplicity vs Layer cluster size in Tracksters;LayerCluster multiplicity in Tracksters;Cluster "
1030  "size (n_{hits})",
1031  nintMplofLCs_,
1032  minMplofLCs_,
1033  maxMplofLCs_,
1034  nintSizeCLsinTSTs_,
1035  minSizeCLsinTSTs_,
1036  maxSizeCLsinTSTs_));
1037 
1038  histograms.h_multiplicity_numberOfEventsHistogram.push_back(ibook.book1D("multiplicity_numberOfEventsHistogram",
1039  "multiplicity numberOfEventsHistogram",
1040  nintMplofLCs_,
1041  minMplofLCs_,
1042  maxMplofLCs_));
1043 
1044  histograms.h_multiplicity_zminus_numberOfEventsHistogram.push_back(
1045  ibook.book1D("multiplicity_zminus_numberOfEventsHistogram",
1046  "multiplicity numberOfEventsHistogram in z-",
1047  nintMplofLCs_,
1048  minMplofLCs_,
1049  maxMplofLCs_));
1050 
1051  histograms.h_multiplicity_zplus_numberOfEventsHistogram.push_back(
1052  ibook.book1D("multiplicity_zplus_numberOfEventsHistogram",
1053  "multiplicity numberOfEventsHistogram in z+",
1054  nintMplofLCs_,
1055  minMplofLCs_,
1056  maxMplofLCs_));
1057 
1058  histograms.h_multiplicityOfLCinTST_vs_layercluster_zminus.push_back(
1059  ibook.book2D("multiplicityOfLCinTST_vs_layercluster_zminus",
1060  "Multiplicity vs Layer number in z-;LayerCluster multiplicity in Tracksters;layer number",
1061  nintMplofLCs_,
1062  minMplofLCs_,
1063  maxMplofLCs_,
1064  layers,
1065  0.,
1066  (float)layers));
1067 
1068  histograms.h_multiplicityOfLCinTST_vs_layercluster_zplus.push_back(
1069  ibook.book2D("multiplicityOfLCinTST_vs_layercluster_zplus",
1070  "Multiplicity vs Layer number in z+;LayerCluster multiplicity in Tracksters;layer number",
1071  nintMplofLCs_,
1072  minMplofLCs_,
1073  maxMplofLCs_,
1074  layers,
1075  0.,
1076  (float)layers));
1077 
1078  histograms.h_multiplicityOfLCinTST_vs_layerclusterenergy.push_back(
1079  ibook.book2D("multiplicityOfLCinTST_vs_layerclusterenergy",
1080  "Multiplicity vs Layer cluster energy;LayerCluster multiplicity in Tracksters;Cluster energy [GeV]",
1081  nintMplofLCs_,
1082  minMplofLCs_,
1083  maxMplofLCs_,
1084  nintClEnepermultiplicity_,
1085  minClEnepermultiplicity_,
1086  maxClEnepermultiplicity_));
1087 
1088  histograms.h_trackster_pt.push_back(
1089  ibook.book1D("trackster_pt", "Pt of the Trackster;Trackster p_{T} [GeV]", nintPt_, minPt_, maxPt_));
1090  histograms.h_trackster_eta.push_back(
1091  ibook.book1D("trackster_eta", "Eta of the Trackster;Trackster #eta", nintEta_, minEta_, maxEta_));
1092  histograms.h_trackster_phi.push_back(
1093  ibook.book1D("trackster_phi", "Phi of the Trackster;Trackster #phi", nintPhi_, minPhi_, maxPhi_));
1094  histograms.h_trackster_energy.push_back(
1095  ibook.book1D("trackster_energy", "Energy of the Trackster;Trackster energy [GeV]", nintEne_, minEne_, maxEne_));
1096  histograms.h_trackster_x.push_back(
1097  ibook.book1D("trackster_x", "X position of the Trackster;Trackster x", nintX_, minX_, maxX_));
1098  histograms.h_trackster_y.push_back(
1099  ibook.book1D("trackster_y", "Y position of the Trackster;Trackster y", nintY_, minY_, maxY_));
1100  histograms.h_trackster_z.push_back(
1101  ibook.book1D("trackster_z", "Z position of the Trackster;Trackster z", nintZ_, minZ_, maxZ_));
1102  histograms.h_trackster_firstlayer.push_back(ibook.book1D(
1103  "trackster_firstlayer", "First layer of the Trackster;Trackster First Layer", 2 * layers, 0., (float)2 * layers));
1104  histograms.h_trackster_lastlayer.push_back(ibook.book1D(
1105  "trackster_lastlayer", "Last layer of the Trackster;Trackster Last Layer", 2 * layers, 0., (float)2 * layers));
1106  histograms.h_trackster_layersnum.push_back(
1107  ibook.book1D("trackster_layersnum",
1108  "Number of layers of the Trackster;Trackster Number of Layers",
1109  2 * layers,
1110  0.,
1111  (float)2 * layers));
1112 }
1113 
1114 void HGVHistoProducerAlgo::bookTracksterSTSHistos(DQMStore::IBooker& ibook,
1115  Histograms& histograms,
1116  const validationType valType) {
1117  const string ref[] = {"caloparticle", "simtrackster", "simtrackster_fromCP"};
1118  const string refT[] = {"CaloParticle", "SimTrackster", "SimTrackster_fromCP"};
1119  // Must be in sync with labels in PostProcessorHGCAL_cfi.py
1120  const string val[] = {"_Link", "_PR"};
1121  //const string val[] = {"_CP", "_PR", "_Link"};
1122  const string rtos = ";score Reco-to-Sim";
1123  const string stor = ";score Sim-to-Reco";
1124  const string shREnFr = ";shared Reco energy fraction";
1125  const string shSEnFr = ";shared Sim energy fraction";
1126 
1127  histograms.h_score_trackster2caloparticle[valType].push_back(
1128  ibook.book1D("Score_trackster2" + ref[valType],
1129  "Score of Trackster per " + refT[valType] + rtos,
1130  nintScore_,
1131  minScore_,
1132  maxScore_));
1133  histograms.h_score_trackster2bestCaloparticle[valType].push_back(
1134  ibook.book1D("ScoreFake_trackster2" + ref[valType],
1135  "Score of Trackster per best " + refT[valType] + rtos,
1136  nintScore_,
1137  minScore_,
1138  maxScore_));
1139  histograms.h_score_trackster2bestCaloparticle2[valType].push_back(
1140  ibook.book1D("ScoreMerge_trackster2" + ref[valType],
1141  "Score of Trackster per 2^{nd} best " + refT[valType] + rtos,
1142  nintScore_,
1143  minScore_,
1144  maxScore_));
1145  histograms.h_score_caloparticle2trackster[valType].push_back(
1146  ibook.book1D("Score_" + ref[valType] + "2trackster",
1147  "Score of " + refT[valType] + " per Trackster" + stor,
1148  nintScore_,
1149  minScore_,
1150  maxScore_));
1151  histograms.h_scorePur_caloparticle2trackster[valType].push_back(
1152  ibook.book1D("ScorePur_" + ref[valType] + "2trackster",
1153  "Score of " + refT[valType] + " per best Trackster" + stor,
1154  nintScore_,
1155  minScore_,
1156  maxScore_));
1157  histograms.h_scoreDupl_caloparticle2trackster[valType].push_back(
1158  ibook.book1D("ScoreDupl_" + ref[valType] + "2trackster",
1159  "Score of " + refT[valType] + " per 2^{nd} best Trackster" + stor,
1160  nintScore_,
1161  minScore_,
1162  maxScore_));
1163  histograms.h_energy_vs_score_trackster2caloparticle[valType].push_back(
1164  ibook.book2D("Energy_vs_Score_trackster2" + refT[valType],
1165  "Energy vs Score of Trackster per " + refT[valType] + rtos + shREnFr,
1166  nintScore_,
1167  minScore_,
1168  maxScore_,
1169  nintSharedEneFrac_,
1170  minTSTSharedEneFrac_,
1171  maxTSTSharedEneFrac_));
1172  histograms.h_energy_vs_score_trackster2bestCaloparticle[valType].push_back(
1173  ibook.book2D("Energy_vs_Score_trackster2best" + refT[valType],
1174  "Energy vs Score of Trackster per best " + refT[valType] + rtos + shREnFr,
1175  nintScore_,
1176  minScore_,
1177  maxScore_,
1178  nintSharedEneFrac_,
1179  minTSTSharedEneFrac_,
1180  maxTSTSharedEneFrac_));
1181  histograms.h_energy_vs_score_trackster2bestCaloparticle2[valType].push_back(
1182  ibook.book2D("Energy_vs_Score_trackster2secBest" + refT[valType],
1183  "Energy vs Score of Trackster per 2^{nd} best " + refT[valType] + rtos + shREnFr,
1184  nintScore_,
1185  minScore_,
1186  maxScore_,
1187  nintSharedEneFrac_,
1188  minTSTSharedEneFrac_,
1189  maxTSTSharedEneFrac_));
1190  histograms.h_energy_vs_score_caloparticle2trackster[valType].push_back(
1191  ibook.book2D("Energy_vs_Score_" + ref[valType] + "2Trackster",
1192  "Energy vs Score of " + refT[valType] + " per Trackster" + stor + shSEnFr,
1193  nintScore_,
1194  minScore_,
1195  maxScore_,
1196  nintSharedEneFrac_,
1197  minTSTSharedEneFrac_,
1198  maxTSTSharedEneFrac_));
1199  histograms.h_energy_vs_score_caloparticle2bestTrackster[valType].push_back(
1200  ibook.book2D("Energy_vs_Score_" + ref[valType] + "2bestTrackster",
1201  "Energy vs Score of " + refT[valType] + " per best Trackster" + stor + shSEnFr,
1202  nintScore_,
1203  minScore_,
1204  maxScore_,
1205  nintSharedEneFrac_,
1206  minTSTSharedEneFrac_,
1207  maxTSTSharedEneFrac_));
1208  histograms.h_energy_vs_score_caloparticle2bestTrackster2[valType].push_back(
1209  ibook.book2D("Energy_vs_Score_" + ref[valType] + "2secBestTrackster",
1210  "Energy vs Score of " + refT[valType] + " per 2^{nd} best Trackster" + stor + shSEnFr,
1211  nintScore_,
1212  minScore_,
1213  maxScore_,
1214  nintSharedEneFrac_,
1215  minTSTSharedEneFrac_,
1216  maxTSTSharedEneFrac_));
1217 
1218  // Back to all Tracksters
1219  // eta
1220  histograms.h_num_trackster_eta[valType].push_back(ibook.book1D(
1221  "Num_Trackster_Eta" + val[valType], "Num Trackster Eta per Trackster;#eta", nintEta_, minEta_, maxEta_));
1222  histograms.h_numMerge_trackster_eta[valType].push_back(ibook.book1D("NumMerge_Trackster_Eta" + val[valType],
1223  "Num Merge Trackster Eta per Trackster;#eta",
1224  nintEta_,
1225  minEta_,
1226  maxEta_));
1227  histograms.h_denom_trackster_eta[valType].push_back(ibook.book1D(
1228  "Denom_Trackster_Eta" + val[valType], "Denom Trackster Eta per Trackster;#eta", nintEta_, minEta_, maxEta_));
1229  // phi
1230  histograms.h_num_trackster_phi[valType].push_back(ibook.book1D(
1231  "Num_Trackster_Phi" + val[valType], "Num Trackster Phi per Trackster;#phi", nintPhi_, minPhi_, maxPhi_));
1232  histograms.h_numMerge_trackster_phi[valType].push_back(ibook.book1D("NumMerge_Trackster_Phi" + val[valType],
1233  "Num Merge Trackster Phi per Trackster;#phi",
1234  nintPhi_,
1235  minPhi_,
1236  maxPhi_));
1237  histograms.h_denom_trackster_phi[valType].push_back(ibook.book1D(
1238  "Denom_Trackster_Phi" + val[valType], "Denom Trackster Phi per Trackster;#phi", nintPhi_, minPhi_, maxPhi_));
1239  // energy
1240  histograms.h_num_trackster_en[valType].push_back(ibook.book1D("Num_Trackster_Energy" + val[valType],
1241  "Num Trackster Energy per Trackster;energy [GeV]",
1242  nintEne_,
1243  minEne_,
1244  maxEne_));
1245  histograms.h_numMerge_trackster_en[valType].push_back(
1246  ibook.book1D("NumMerge_Trackster_Energy" + val[valType],
1247  "Num Merge Trackster Energy per Trackster;energy [GeV]",
1248  nintEne_,
1249  minEne_,
1250  maxEne_));
1251  histograms.h_denom_trackster_en[valType].push_back(ibook.book1D("Denom_Trackster_Energy" + val[valType],
1252  "Denom Trackster Energy per Trackster;energy [GeV]",
1253  nintEne_,
1254  minEne_,
1255  maxEne_));
1256  // pT
1257  histograms.h_num_trackster_pt[valType].push_back(ibook.book1D(
1258  "Num_Trackster_Pt" + val[valType], "Num Trackster p_{T} per Trackster;p_{T} [GeV]", nintPt_, minPt_, maxPt_));
1259  histograms.h_numMerge_trackster_pt[valType].push_back(
1260  ibook.book1D("NumMerge_Trackster_Pt" + val[valType],
1261  "Num Merge Trackster p_{T} per Trackster;p_{T} [GeV]",
1262  nintPt_,
1263  minPt_,
1264  maxPt_));
1265  histograms.h_denom_trackster_pt[valType].push_back(ibook.book1D(
1266  "Denom_Trackster_Pt" + val[valType], "Denom Trackster p_{T} per Trackster;p_{T} [GeV]", nintPt_, minPt_, maxPt_));
1267 
1268  histograms.h_sharedenergy_trackster2caloparticle[valType].push_back(
1269  ibook.book1D("SharedEnergy_trackster2" + ref[valType],
1270  "Shared Energy of Trackster per " + refT[valType] + shREnFr,
1271  nintSharedEneFrac_,
1272  minTSTSharedEneFrac_,
1273  maxTSTSharedEneFrac_));
1274  histograms.h_sharedenergy_trackster2bestCaloparticle[valType].push_back(
1275  ibook.book1D("SharedEnergy_trackster2" + ref[valType] + "_assoc",
1276  "Shared Energy of Trackster per best " + refT[valType] + shREnFr,
1277  nintSharedEneFrac_,
1278  minTSTSharedEneFrac_,
1279  maxTSTSharedEneFrac_));
1280  histograms.h_sharedenergy_trackster2bestCaloparticle_vs_eta[valType].push_back(
1281  ibook.bookProfile("SharedEnergy_trackster2" + ref[valType] + "_assoc_vs_eta",
1282  "Shared Energy of Trackster vs #eta per best " + refT[valType] + ";Trackster #eta" + shREnFr,
1283  nintEta_,
1284  minEta_,
1285  maxEta_,
1286  minTSTSharedEneFrac_,
1287  maxTSTSharedEneFrac_));
1288  histograms.h_sharedenergy_trackster2bestCaloparticle_vs_phi[valType].push_back(
1289  ibook.bookProfile("SharedEnergy_trackster2" + ref[valType] + "_assoc_vs_phi",
1290  "Shared Energy of Trackster vs #phi per best " + refT[valType] + ";Trackster #phi" + shREnFr,
1291  nintPhi_,
1292  minPhi_,
1293  maxPhi_,
1294  minTSTSharedEneFrac_,
1295  maxTSTSharedEneFrac_));
1296  histograms.h_sharedenergy_trackster2bestCaloparticle2[valType].push_back(
1297  ibook.book1D("SharedEnergy_trackster2" + ref[valType] + "_assoc2",
1298  "Shared Energy of Trackster per 2^{nd} best " + refT[valType] + shREnFr,
1299  nintSharedEneFrac_,
1300  minTSTSharedEneFrac_,
1301  maxTSTSharedEneFrac_));
1302 
1303  histograms.h_sharedenergy_caloparticle2trackster[valType].push_back(
1304  ibook.book1D("SharedEnergy_" + ref[valType] + "2trackster",
1305  "Shared Energy of " + refT[valType] + " per Trackster" + shSEnFr,
1306  nintSharedEneFrac_,
1307  minTSTSharedEneFrac_,
1308  maxTSTSharedEneFrac_));
1309  histograms.h_sharedenergy_caloparticle2trackster_assoc[valType].push_back(
1310  ibook.book1D("SharedEnergy_" + ref[valType] + "2trackster_assoc",
1311  "Shared Energy of " + refT[valType] + " per best Trackster" + shSEnFr,
1312  nintSharedEneFrac_,
1313  minTSTSharedEneFrac_,
1314  maxTSTSharedEneFrac_));
1315  histograms.h_sharedenergy_caloparticle2trackster_assoc_vs_eta[valType].push_back(ibook.bookProfile(
1316  "SharedEnergy_" + ref[valType] + "2trackster_assoc_vs_eta",
1317  "Shared Energy of " + refT[valType] + " vs #eta per best Trackster;" + refT[valType] + " #eta" + shSEnFr,
1318  nintEta_,
1319  minEta_,
1320  maxEta_,
1321  minTSTSharedEneFrac_,
1322  maxTSTSharedEneFrac_));
1323  histograms.h_sharedenergy_caloparticle2trackster_assoc_vs_phi[valType].push_back(ibook.bookProfile(
1324  "SharedEnergy_" + ref[valType] + "2trackster_assoc_vs_phi",
1325  "Shared Energy of " + refT[valType] + " vs #phi per best Trackster;" + refT[valType] + " #phi" + shSEnFr,
1326  nintPhi_,
1327  minPhi_,
1328  maxPhi_,
1329  minTSTSharedEneFrac_,
1330  maxTSTSharedEneFrac_));
1331  histograms.h_sharedenergy_caloparticle2trackster_assoc2[valType].push_back(
1332  ibook.book1D("SharedEnergy_" + ref[valType] + "2trackster_assoc2",
1333  "Shared Energy of " + refT[valType] + " per 2^{nd} best Trackster;" + shSEnFr,
1334  nintSharedEneFrac_,
1335  minTSTSharedEneFrac_,
1336  maxTSTSharedEneFrac_));
1337 
1338  // eta
1339  histograms.h_numEff_caloparticle_eta[valType].push_back(
1340  ibook.book1D("NumEff_" + refT[valType] + "_Eta",
1341  "Num Efficiency " + refT[valType] + " Eta per Trackster;#eta",
1342  nintEta_,
1343  minEta_,
1344  maxEta_));
1345  histograms.h_num_caloparticle_eta[valType].push_back(
1346  ibook.book1D("Num_" + refT[valType] + "_Eta",
1347  "Num Purity " + refT[valType] + " Eta per Trackster;#eta",
1348  nintEta_,
1349  minEta_,
1350  maxEta_));
1351  histograms.h_numDup_trackster_eta[valType].push_back(ibook.book1D(
1352  "NumDup_Trackster_Eta" + val[valType], "Num Duplicate Trackster vs Eta;#eta", nintEta_, minEta_, maxEta_));
1353  histograms.h_denom_caloparticle_eta[valType].push_back(
1354  ibook.book1D("Denom_" + refT[valType] + "_Eta",
1355  "Denom " + refT[valType] + " Eta per Trackster;#eta",
1356  nintEta_,
1357  minEta_,
1358  maxEta_));
1359  // phi
1360  histograms.h_numEff_caloparticle_phi[valType].push_back(
1361  ibook.book1D("NumEff_" + refT[valType] + "_Phi",
1362  "Num Efficiency " + refT[valType] + " Phi per Trackster;#phi",
1363  nintPhi_,
1364  minPhi_,
1365  maxPhi_));
1366  histograms.h_num_caloparticle_phi[valType].push_back(
1367  ibook.book1D("Num_" + refT[valType] + "_Phi",
1368  "Num Purity " + refT[valType] + " Phi per Trackster;#phi",
1369  nintPhi_,
1370  minPhi_,
1371  maxPhi_));
1372  histograms.h_numDup_trackster_phi[valType].push_back(ibook.book1D(
1373  "NumDup_Trackster_Phi" + val[valType], "Num Duplicate Trackster vs Phi;#phi", nintPhi_, minPhi_, maxPhi_));
1374  histograms.h_denom_caloparticle_phi[valType].push_back(
1375  ibook.book1D("Denom_" + refT[valType] + "_Phi",
1376  "Denom " + refT[valType] + " Phi per Trackster;#phi",
1377  nintPhi_,
1378  minPhi_,
1379  maxPhi_));
1380  // energy
1381  histograms.h_numEff_caloparticle_en[valType].push_back(
1382  ibook.book1D("NumEff_" + refT[valType] + "_Energy",
1383  "Num Efficiency " + refT[valType] + " Energy per Trackster;energy [GeV]",
1384  nintEne_,
1385  minEne_,
1386  maxEne_));
1387  histograms.h_num_caloparticle_en[valType].push_back(
1388  ibook.book1D("Num_" + refT[valType] + "_Energy",
1389  "Num Purity " + refT[valType] + " Energy per Trackster;energy [GeV]",
1390  nintEne_,
1391  minEne_,
1392  maxEne_));
1393  histograms.h_numDup_trackster_en[valType].push_back(ibook.book1D("NumDup_Trackster_Energy" + val[valType],
1394  "Num Duplicate Trackster vs Energy;energy [GeV]",
1395  nintEne_,
1396  minEne_,
1397  maxEne_));
1398  histograms.h_denom_caloparticle_en[valType].push_back(
1399  ibook.book1D("Denom_" + refT[valType] + "_Energy",
1400  "Denom " + refT[valType] + " Energy per Trackster;energy [GeV]",
1401  nintEne_,
1402  minEne_,
1403  maxEne_));
1404  // pT
1405  histograms.h_numEff_caloparticle_pt[valType].push_back(
1406  ibook.book1D("NumEff_" + refT[valType] + "_Pt",
1407  "Num Efficiency " + refT[valType] + " p_{T} per Trackster;p_{T} [GeV]",
1408  nintPt_,
1409  minPt_,
1410  maxPt_));
1411  histograms.h_num_caloparticle_pt[valType].push_back(
1412  ibook.book1D("Num_" + refT[valType] + "_Pt",
1413  "Num Purity " + refT[valType] + " p_{T} per Trackster;p_{T} [GeV]",
1414  nintPt_,
1415  minPt_,
1416  maxPt_));
1417  histograms.h_numDup_trackster_pt[valType].push_back(ibook.book1D(
1418  "NumDup_Trackster_Pt" + val[valType], "Num Duplicate Trackster vs p_{T};p_{T} [GeV]", nintPt_, minPt_, maxPt_));
1419  histograms.h_denom_caloparticle_pt[valType].push_back(
1420  ibook.book1D("Denom_" + refT[valType] + "_Pt",
1421  "Denom " + refT[valType] + " p_{T} per Trackster;p_{T} [GeV]",
1422  nintPt_,
1423  minPt_,
1424  maxPt_));
1425 }
1426 
1427 void HGVHistoProducerAlgo::fill_info_histos(const Histograms& histograms, unsigned int layers) const {
1428  // Save some info straight from geometry to avoid mistakes from updates
1429  //----------- TODO ----------------------------------------------------------
1430  // For now values returned for 'lastLayerFHzp': '104', 'lastLayerFHzm': '52' are not the one expected.
1431  // Will come back to this when there will be info in CMSSW to put in DQM file.
1432  histograms.lastLayerEEzm->Fill(recHitTools_->lastLayerEE());
1433  histograms.lastLayerFHzm->Fill(recHitTools_->lastLayerFH());
1434  histograms.maxlayerzm->Fill(layers);
1435  histograms.lastLayerEEzp->Fill(recHitTools_->lastLayerEE() + layers);
1436  histograms.lastLayerFHzp->Fill(recHitTools_->lastLayerFH() + layers);
1437  histograms.maxlayerzp->Fill(layers + layers);
1438 }
1439 
1440 void HGVHistoProducerAlgo::fill_caloparticle_histos(const Histograms& histograms,
1441  int pdgid,
1442  const CaloParticle& caloParticle,
1443  std::vector<SimVertex> const& simVertices,
1444  unsigned int layers,
1445  std::unordered_map<DetId, const HGCRecHit*> const& hitMap) const {
1446  const auto eta = getEta(caloParticle.eta());
1447  if (histograms.h_caloparticle_eta.count(pdgid)) {
1448  histograms.h_caloparticle_eta.at(pdgid)->Fill(eta);
1449  }
1450  if (histograms.h_caloparticle_eta_Zorigin.count(pdgid)) {
1451  histograms.h_caloparticle_eta_Zorigin.at(pdgid)->Fill(
1452  simVertices.at(caloParticle.g4Tracks()[0].vertIndex()).position().z(), eta);
1453  }
1454 
1455  if (histograms.h_caloparticle_energy.count(pdgid)) {
1456  histograms.h_caloparticle_energy.at(pdgid)->Fill(caloParticle.energy());
1457  }
1458  if (histograms.h_caloparticle_pt.count(pdgid)) {
1459  histograms.h_caloparticle_pt.at(pdgid)->Fill(caloParticle.pt());
1460  }
1461  if (histograms.h_caloparticle_phi.count(pdgid)) {
1462  histograms.h_caloparticle_phi.at(pdgid)->Fill(caloParticle.phi());
1463  }
1464 
1465  if (histograms.h_caloparticle_nSimClusters.count(pdgid)) {
1466  histograms.h_caloparticle_nSimClusters.at(pdgid)->Fill(caloParticle.simClusters().size());
1467 
1468  int simHits = 0;
1469  int minLayerId = 999;
1470  int maxLayerId = 0;
1471 
1472  int simHits_matched = 0;
1473  int minLayerId_matched = 999;
1474  int maxLayerId_matched = 0;
1475 
1476  float energy = 0.;
1477  std::map<int, double> totenergy_layer;
1478 
1479  float hitEnergyWeight_invSum = 0;
1480  std::vector<std::pair<DetId, float>> haf_cp;
1481  for (const auto& sc : caloParticle.simClusters()) {
1482  LogDebug("HGCalValidator") << " This sim cluster has " << sc->hits_and_fractions().size() << " simHits and "
1483  << sc->energy() << " energy. " << std::endl;
1484  simHits += sc->hits_and_fractions().size();
1485  for (auto const& h_and_f : sc->hits_and_fractions()) {
1486  const auto hitDetId = h_and_f.first;
1487  const int layerId =
1488  recHitTools_->getLayerWithOffset(hitDetId) + layers * ((recHitTools_->zside(hitDetId) + 1) >> 1) - 1;
1489  // set to 0 if matched RecHit not found
1490  int layerId_matched_min = 999;
1491  int layerId_matched_max = 0;
1492  std::unordered_map<DetId, const HGCRecHit*>::const_iterator itcheck = hitMap.find(hitDetId);
1493  if (itcheck != hitMap.end()) {
1494  layerId_matched_min = layerId;
1495  layerId_matched_max = layerId;
1496  simHits_matched++;
1497 
1498  const auto hitEn = itcheck->second->energy();
1499  hitEnergyWeight_invSum += pow(hitEn, 2);
1500  const auto hitFr = h_and_f.second;
1501  const auto hitEnFr = hitEn * hitFr;
1502  energy += hitEnFr;
1503  histograms.h_caloparticle_nHits_matched_energy.at(pdgid)->Fill(hitEnFr);
1504  histograms.h_caloparticle_nHits_matched_energy_layer.at(pdgid)->Fill(layerId, hitEnFr);
1505 
1506  if (totenergy_layer.find(layerId) != totenergy_layer.end()) {
1507  totenergy_layer[layerId] = totenergy_layer.at(layerId) + hitEn;
1508  } else {
1509  totenergy_layer.emplace(layerId, hitEn);
1510  }
1511  if (caloParticle.simClusters().size() == 1)
1512  histograms.h_caloparticle_nHits_matched_energy_layer_1SimCl.at(pdgid)->Fill(layerId, hitEnFr);
1513 
1514  auto found = std::find_if(std::begin(haf_cp),
1515  std::end(haf_cp),
1516  [&hitDetId](const std::pair<DetId, float>& v) { return v.first == hitDetId; });
1517  if (found != haf_cp.end())
1518  found->second += hitFr;
1519  else
1520  haf_cp.emplace_back(hitDetId, hitFr);
1521 
1522  } else {
1523  LogDebug("HGCalValidator") << " matched to RecHit NOT found !" << std::endl;
1524  }
1525 
1526  minLayerId = std::min(minLayerId, layerId);
1527  maxLayerId = std::max(maxLayerId, layerId);
1528  minLayerId_matched = std::min(minLayerId_matched, layerId_matched_min);
1529  maxLayerId_matched = std::max(maxLayerId_matched, layerId_matched_max);
1530  }
1531  LogDebug("HGCalValidator") << std::endl;
1532  } // End loop over SimClusters of CaloParticle
1533  if (hitEnergyWeight_invSum)
1534  hitEnergyWeight_invSum = 1 / hitEnergyWeight_invSum;
1535 
1536  histograms.h_caloparticle_firstlayer.at(pdgid)->Fill(minLayerId);
1537  histograms.h_caloparticle_lastlayer.at(pdgid)->Fill(maxLayerId);
1538  histograms.h_caloparticle_layersnum.at(pdgid)->Fill(int(maxLayerId - minLayerId));
1539 
1540  histograms.h_caloparticle_firstlayer_matchedtoRecHit.at(pdgid)->Fill(minLayerId_matched);
1541  histograms.h_caloparticle_lastlayer_matchedtoRecHit.at(pdgid)->Fill(maxLayerId_matched);
1542  histograms.h_caloparticle_layersnum_matchedtoRecHit.at(pdgid)->Fill(int(maxLayerId_matched - minLayerId_matched));
1543 
1544  histograms.h_caloparticle_nHitsInSimClusters.at(pdgid)->Fill((float)simHits);
1545  histograms.h_caloparticle_nHitsInSimClusters_matchedtoRecHit.at(pdgid)->Fill((float)simHits_matched);
1546  histograms.h_caloparticle_selfenergy.at(pdgid)->Fill((float)energy);
1547  histograms.h_caloparticle_energyDifference.at(pdgid)->Fill((float)1. - energy / caloParticle.energy());
1548 
1549  //Calculate sum energy per-layer
1550  auto i = totenergy_layer.begin();
1551  double sum_energy = 0.0;
1552  while (i != totenergy_layer.end()) {
1553  sum_energy += i->second;
1554  histograms.h_caloparticle_sum_energy_layer.at(pdgid)->Fill(i->first, sum_energy / caloParticle.energy() * 100.);
1555  i++;
1556  }
1557 
1558  for (auto const& haf : haf_cp) {
1559  const auto hitEn = hitMap.find(haf.first)->second->energy();
1560  const auto weight = pow(hitEn, 2);
1561  histograms.h_caloparticle_fractions.at(pdgid)->Fill(haf.second, weight * hitEnergyWeight_invSum);
1562  histograms.h_caloparticle_fractions_weight.at(pdgid)->Fill(haf.second, weight * hitEnergyWeight_invSum, weight);
1563  }
1564  }
1565 }
1566 
1567 void HGVHistoProducerAlgo::HGVHistoProducerAlgo::fill_simCluster_histos(const Histograms& histograms,
1568  std::vector<SimCluster> const& simClusters,
1569  unsigned int layers,
1570  std::vector<int> thicknesses) const {
1571  //Each event to be treated as two events: an event in +ve endcap,
1572  //plus another event in -ve endcap. In this spirit there will be
1573  //a layer variable (layerid) that maps the layers in :
1574  //-z: 0->49
1575  //+z: 50->99
1576 
1577  //To keep track of total num of simClusters per layer
1578  //tnscpl[layerid]
1579  std::vector<int> tnscpl(1000, 0); //tnscpl.clear(); tnscpl.reserve(1000);
1580 
1581  //To keep track of the total num of clusters per thickness in plus and in minus endcaps
1582  std::map<std::string, int> tnscpthplus;
1583  tnscpthplus.clear();
1584  std::map<std::string, int> tnscpthminus;
1585  tnscpthminus.clear();
1586  //At the beginning of the event all layers should be initialized to zero total clusters per thickness
1587  for (std::vector<int>::iterator it = thicknesses.begin(); it != thicknesses.end(); ++it) {
1588  tnscpthplus.insert(std::pair<std::string, int>(std::to_string(*it), 0));
1589  tnscpthminus.insert(std::pair<std::string, int>(std::to_string(*it), 0));
1590  }
1591  //To keep track of the total num of simClusters with mixed thickness hits per event
1592  tnscpthplus.insert(std::pair<std::string, int>("mixed", 0));
1593  tnscpthminus.insert(std::pair<std::string, int>("mixed", 0));
1594 
1595  //loop through simClusters
1596  for (const auto& sc : simClusters) {
1597  //Auxillary variables to count the number of different kind of hits in each simCluster
1598  int nthhits120p = 0;
1599  int nthhits200p = 0;
1600  int nthhits300p = 0;
1601  int nthhitsscintp = 0;
1602  int nthhits120m = 0;
1603  int nthhits200m = 0;
1604  int nthhits300m = 0;
1605  int nthhitsscintm = 0;
1606  //For the hits thickness of the layer cluster.
1607  double thickness = 0.;
1608  //To keep track if we added the simCluster in a specific layer
1609  std::vector<int> occurenceSCinlayer(1000, 0); //[layerid][0 if not added]
1610 
1611  //loop through hits of the simCluster
1612  for (const auto& hAndF : sc.hits_and_fractions()) {
1613  const DetId sh_detid = hAndF.first;
1614 
1615  //The layer the cluster belongs to. As mentioned in the mapping above, it takes into account -z and +z.
1616  int layerid =
1617  recHitTools_->getLayerWithOffset(sh_detid) + layers * ((recHitTools_->zside(sh_detid) + 1) >> 1) - 1;
1618  //zside that the current cluster belongs to.
1619  int zside = recHitTools_->zside(sh_detid);
1620 
1621  //add the simCluster to the relevant layer. A SimCluster may give contribution to several layers.
1622  if (occurenceSCinlayer[layerid] == 0) {
1623  tnscpl[layerid]++;
1624  }
1625  occurenceSCinlayer[layerid]++;
1626 
1627  if (sh_detid.det() == DetId::Forward || sh_detid.det() == DetId::HGCalEE || sh_detid.det() == DetId::HGCalHSi) {
1628  thickness = recHitTools_->getSiThickness(sh_detid);
1629  } else if (sh_detid.det() == DetId::HGCalHSc) {
1630  thickness = -1;
1631  } else {
1632  LogDebug("HGCalValidator") << "These are HGCal simClusters, you shouldn't be here !!! " << layerid << "\n";
1633  continue;
1634  }
1635 
1636  if ((thickness == 120.) && (zside > 0.)) {
1637  nthhits120p++;
1638  } else if ((thickness == 120.) && (zside < 0.)) {
1639  nthhits120m++;
1640  } else if ((thickness == 200.) && (zside > 0.)) {
1641  nthhits200p++;
1642  } else if ((thickness == 200.) && (zside < 0.)) {
1643  nthhits200m++;
1644  } else if ((thickness == 300.) && (zside > 0.)) {
1645  nthhits300p++;
1646  } else if ((thickness == 300.) && (zside < 0.)) {
1647  nthhits300m++;
1648  } else if ((thickness == -1) && (zside > 0.)) {
1649  nthhitsscintp++;
1650  } else if ((thickness == -1) && (zside < 0.)) {
1651  nthhitsscintm++;
1652  } else { //assert(0);
1653  LogDebug("HGCalValidator")
1654  << " You are running a geometry that contains thicknesses different than the normal ones. "
1655  << "\n";
1656  }
1657 
1658  } //end of loop through hits
1659 
1660  //Check for simultaneously having hits of different kind. Checking at least two combinations is sufficient.
1661  if ((nthhits120p != 0 && nthhits200p != 0) || (nthhits120p != 0 && nthhits300p != 0) ||
1662  (nthhits120p != 0 && nthhitsscintp != 0) || (nthhits200p != 0 && nthhits300p != 0) ||
1663  (nthhits200p != 0 && nthhitsscintp != 0) || (nthhits300p != 0 && nthhitsscintp != 0)) {
1664  tnscpthplus["mixed"]++;
1665  } else if ((nthhits120p != 0 || nthhits200p != 0 || nthhits300p != 0 || nthhitsscintp != 0)) {
1666  //This is a cluster with hits of one kind
1667  tnscpthplus[std::to_string((int)thickness)]++;
1668  }
1669  if ((nthhits120m != 0 && nthhits200m != 0) || (nthhits120m != 0 && nthhits300m != 0) ||
1670  (nthhits120m != 0 && nthhitsscintm != 0) || (nthhits200m != 0 && nthhits300m != 0) ||
1671  (nthhits200m != 0 && nthhitsscintm != 0) || (nthhits300m != 0 && nthhitsscintm != 0)) {
1672  tnscpthminus["mixed"]++;
1673  } else if ((nthhits120m != 0 || nthhits200m != 0 || nthhits300m != 0 || nthhitsscintm != 0)) {
1674  //This is a cluster with hits of one kind
1675  tnscpthminus[std::to_string((int)thickness)]++;
1676  }
1677 
1678  } //end of loop through SimClusters of the event
1679 
1680  //Per layer : Loop 0->99
1681  for (unsigned ilayer = 0; ilayer < layers * 2; ++ilayer)
1682  if (histograms.h_simclusternum_perlayer.count(ilayer))
1683  histograms.h_simclusternum_perlayer.at(ilayer)->Fill(tnscpl[ilayer]);
1684 
1685  //Per thickness
1686  for (std::vector<int>::iterator it = thicknesses.begin(); it != thicknesses.end(); ++it) {
1687  if (histograms.h_simclusternum_perthick.count(*it)) {
1688  histograms.h_simclusternum_perthick.at(*it)->Fill(tnscpthplus[std::to_string(*it)]);
1689  histograms.h_simclusternum_perthick.at(*it)->Fill(tnscpthminus[std::to_string(*it)]);
1690  }
1691  }
1692  //Mixed thickness clusters
1693  histograms.h_mixedhitssimcluster_zplus->Fill(tnscpthplus["mixed"]);
1694  histograms.h_mixedhitssimcluster_zminus->Fill(tnscpthminus["mixed"]);
1695 }
1696 
1697 void HGVHistoProducerAlgo::HGVHistoProducerAlgo::fill_simClusterAssociation_histos(
1698  const Histograms& histograms,
1699  const int count,
1700  edm::Handle<reco::CaloClusterCollection> clusterHandle,
1701  const reco::CaloClusterCollection& clusters,
1702  edm::Handle<std::vector<SimCluster>> simClusterHandle,
1703  std::vector<SimCluster> const& simClusters,
1704  std::vector<size_t> const& sCIndices,
1705  const std::vector<float>& mask,
1706  std::unordered_map<DetId, const HGCRecHit*> const& hitMap,
1707  unsigned int layers,
1708  const hgcal::RecoToSimCollectionWithSimClusters& scsInLayerClusterMap,
1709  const hgcal::SimToRecoCollectionWithSimClusters& lcsInSimClusterMap) const {
1710  //Each event to be treated as two events: an event in +ve endcap,
1711  //plus another event in -ve endcap. In this spirit there will be
1712  //a layer variable (layerid) that maps the layers in :
1713  //-z: 0->49
1714  //+z: 50->99
1715 
1716  //Will add some general plots on the specific mask in the future.
1717 
1718  layerClusters_to_SimClusters(histograms,
1719  count,
1720  clusterHandle,
1721  clusters,
1722  simClusterHandle,
1723  simClusters,
1724  sCIndices,
1725  mask,
1726  hitMap,
1727  layers,
1728  scsInLayerClusterMap,
1729  lcsInSimClusterMap);
1730 }
1731 
1732 void HGVHistoProducerAlgo::fill_cluster_histos(const Histograms& histograms,
1733  const int count,
1734  const reco::CaloCluster& cluster) const {
1735  const auto eta = getEta(cluster.eta());
1736  histograms.h_cluster_eta[count]->Fill(eta);
1737 }
1738 
1739 void HGVHistoProducerAlgo::layerClusters_to_CaloParticles(const Histograms& histograms,
1740  edm::Handle<reco::CaloClusterCollection> clusterHandle,
1741  const reco::CaloClusterCollection& clusters,
1742  edm::Handle<std::vector<CaloParticle>> caloParticleHandle,
1743  std::vector<CaloParticle> const& cP,
1744  std::vector<size_t> const& cPIndices,
1745  std::vector<size_t> const& cPSelectedIndices,
1746  std::unordered_map<DetId, const HGCRecHit*> const& hitMap,
1747  unsigned int layers,
1748  const hgcal::RecoToSimCollection& cpsInLayerClusterMap,
1749  const hgcal::SimToRecoCollection& cPOnLayerMap) const {
1750  const auto nLayerClusters = clusters.size();
1751 
1752  std::unordered_map<DetId, std::vector<HGVHistoProducerAlgo::detIdInfoInCluster>> detIdToCaloParticleId_Map;
1753  std::unordered_map<DetId, std::vector<HGVHistoProducerAlgo::detIdInfoInCluster>> detIdToLayerClusterId_Map;
1754 
1755  // The association has to be done in an all-vs-all fashion.
1756  // For this reason use the full set of CaloParticles, with the only filter on bx
1757  for (const auto& cpId : cPIndices) {
1758  for (const auto& simCluster : cP[cpId].simClusters()) {
1759  for (const auto& it_haf : simCluster->hits_and_fractions()) {
1760  const DetId hitid = (it_haf.first);
1761  if (hitMap.find(hitid) != hitMap.end()) {
1762  if (detIdToCaloParticleId_Map.find(hitid) == detIdToCaloParticleId_Map.end()) {
1763  detIdToCaloParticleId_Map[hitid] = std::vector<HGVHistoProducerAlgo::detIdInfoInCluster>();
1764  detIdToCaloParticleId_Map[hitid].emplace_back(
1765  HGVHistoProducerAlgo::detIdInfoInCluster{cpId, it_haf.second});
1766  } else {
1767  auto findHitIt =
1768  std::find(detIdToCaloParticleId_Map[hitid].begin(),
1769  detIdToCaloParticleId_Map[hitid].end(),
1770  HGVHistoProducerAlgo::detIdInfoInCluster{
1771  cpId, 0.f}); // only the first element is used for the matching (overloaded operator==)
1772  if (findHitIt != detIdToCaloParticleId_Map[hitid].end())
1773  findHitIt->fraction += it_haf.second;
1774  else
1775  detIdToCaloParticleId_Map[hitid].emplace_back(
1776  HGVHistoProducerAlgo::detIdInfoInCluster{cpId, it_haf.second});
1777  }
1778  }
1779  }
1780  }
1781  }
1782 
1783  for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) {
1784  const auto& hits_and_fractions = clusters[lcId].hitsAndFractions();
1785  const auto numberOfHitsInLC = hits_and_fractions.size();
1786 
1787  // This vector will store, for each hit in the Layercluster, the index of
1788  // the CaloParticle that contributed the most, in terms of energy, to it.
1789  // Special values are:
1790  //
1791  // -2 --> the reconstruction fraction of the RecHit is 0 (used in the past to monitor Halo Hits)
1792  // -3 --> same as before with the added condition that no CaloParticle has been linked to this RecHit
1793  // -1 --> the reco fraction is >0, but no CaloParticle has been linked to it
1794  // >=0 --> index of the linked CaloParticle
1795  std::vector<int> hitsToCaloParticleId(numberOfHitsInLC);
1796  const auto firstHitDetId = hits_and_fractions[0].first;
1797  int lcLayerId =
1798  recHitTools_->getLayerWithOffset(firstHitDetId) + layers * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
1799 
1800  // This will store the fraction of the CaloParticle energy shared with the LayerCluster: e_shared/cp_energy
1801  std::unordered_map<unsigned, float> CPEnergyInLC;
1802 
1803  for (unsigned int iHit = 0; iHit < numberOfHitsInLC; iHit++) {
1804  const DetId rh_detid = hits_and_fractions[iHit].first;
1805  const auto rhFraction = hits_and_fractions[iHit].second;
1806 
1807  std::unordered_map<DetId, const HGCRecHit*>::const_iterator itcheck = hitMap.find(rh_detid);
1808  const HGCRecHit* hit = itcheck->second;
1809 
1810  if (detIdToLayerClusterId_Map.find(rh_detid) == detIdToLayerClusterId_Map.end()) {
1811  detIdToLayerClusterId_Map[rh_detid] = std::vector<HGVHistoProducerAlgo::detIdInfoInCluster>();
1812  }
1813  detIdToLayerClusterId_Map[rh_detid].emplace_back(HGVHistoProducerAlgo::detIdInfoInCluster{lcId, rhFraction});
1814 
1815  const auto& hit_find_in_CP = detIdToCaloParticleId_Map.find(rh_detid);
1816 
1817  // if the fraction is zero or the hit does not belong to any calo
1818  // particle, set the caloparticleId for the hit to -1 this will
1819  // contribute to the number of noise hits
1820 
1821  // MR Remove the case in which the fraction is 0, since this could be a
1822  // real hit that has been marked as halo.
1823  if (rhFraction == 0.) {
1824  hitsToCaloParticleId[iHit] = -2;
1825  }
1826  if (hit_find_in_CP == detIdToCaloParticleId_Map.end()) {
1827  hitsToCaloParticleId[iHit] -= 1;
1828  } else {
1829  auto maxCPEnergyInLC = 0.f;
1830  auto maxCPId = -1;
1831  for (auto& h : hit_find_in_CP->second) {
1832  const auto iCP = h.clusterId;
1833  CPEnergyInLC[iCP] += h.fraction * hit->energy();
1834  // Keep track of which CaloParticle contributed the most, in terms
1835  // of energy, to this specific LayerCluster.
1836  if (CPEnergyInLC[iCP] > maxCPEnergyInLC) {
1837  maxCPEnergyInLC = CPEnergyInLC[iCP];
1838  maxCPId = iCP;
1839  }
1840  }
1841  hitsToCaloParticleId[iHit] = maxCPId;
1842  }
1843  histograms.h_cellAssociation_perlayer.at(lcLayerId)->Fill(
1844  hitsToCaloParticleId[iHit] > 0. ? 0. : hitsToCaloParticleId[iHit]);
1845  } // End loop over hits on a LayerCluster
1846 
1847  } // End of loop over LayerClusters
1848 
1849  // Fill the plots to compute the different metrics linked to
1850  // reco-level, namely fake-rate an merge-rate. In this loop should *not*
1851  // restrict only to the selected caloParaticles.
1852  for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) {
1853  const auto firstHitDetId = (clusters[lcId].hitsAndFractions())[0].first;
1854  const int lcLayerId =
1855  recHitTools_->getLayerWithOffset(firstHitDetId) + layers * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
1856  histograms.h_denom_layercl_eta_perlayer.at(lcLayerId)->Fill(clusters[lcId].eta());
1857  histograms.h_denom_layercl_phi_perlayer.at(lcLayerId)->Fill(clusters[lcId].phi());
1858  //
1859  const edm::Ref<reco::CaloClusterCollection> lcRef(clusterHandle, lcId);
1860  const auto& cpsIt = cpsInLayerClusterMap.find(lcRef);
1861  if (cpsIt == cpsInLayerClusterMap.end())
1862  continue;
1863 
1864  const auto lc_en = clusters[lcId].energy();
1865  const auto& cps = cpsIt->val;
1866  if (lc_en == 0. && !cps.empty()) {
1867  for (const auto& cpPair : cps)
1868  histograms.h_score_layercl2caloparticle_perlayer.at(lcLayerId)->Fill(cpPair.second);
1869  continue;
1870  }
1871  for (const auto& cpPair : cps) {
1872  LogDebug("HGCalValidator") << "layerCluster Id: \t" << lcId << "\t CP id: \t" << cpPair.first.index()
1873  << "\t score \t" << cpPair.second << std::endl;
1874  histograms.h_score_layercl2caloparticle_perlayer.at(lcLayerId)->Fill(cpPair.second);
1875  auto const& cp_linked =
1876  std::find_if(std::begin(cPOnLayerMap[cpPair.first]),
1877  std::end(cPOnLayerMap[cpPair.first]),
1878  [&lcRef](const std::pair<edm::Ref<reco::CaloClusterCollection>, std::pair<float, float>>& p) {
1879  return p.first == lcRef;
1880  });
1881  if (cp_linked ==
1882  cPOnLayerMap[cpPair.first].end()) // This should never happen by construction of the association maps
1883  continue;
1884  histograms.h_sharedenergy_layercl2caloparticle_perlayer.at(lcLayerId)->Fill(cp_linked->second.first / lc_en,
1885  lc_en);
1886  histograms.h_energy_vs_score_layercl2caloparticle_perlayer.at(lcLayerId)->Fill(cpPair.second,
1887  cp_linked->second.first / lc_en);
1888  }
1889  const auto assoc =
1890  std::count_if(std::begin(cps), std::end(cps), [](const auto& obj) { return obj.second < ScoreCutLCtoCP_; });
1891  if (assoc) {
1892  histograms.h_num_layercl_eta_perlayer.at(lcLayerId)->Fill(clusters[lcId].eta());
1893  histograms.h_num_layercl_phi_perlayer.at(lcLayerId)->Fill(clusters[lcId].phi());
1894  if (assoc > 1) {
1895  histograms.h_numMerge_layercl_eta_perlayer.at(lcLayerId)->Fill(clusters[lcId].eta());
1896  histograms.h_numMerge_layercl_phi_perlayer.at(lcLayerId)->Fill(clusters[lcId].phi());
1897  }
1898  const auto& best = std::min_element(
1899  std::begin(cps), std::end(cps), [](const auto& obj1, const auto& obj2) { return obj1.second < obj2.second; });
1900  const auto& best_cp_linked =
1901  std::find_if(std::begin(cPOnLayerMap[best->first]),
1902  std::end(cPOnLayerMap[best->first]),
1903  [&lcRef](const std::pair<edm::Ref<reco::CaloClusterCollection>, std::pair<float, float>>& p) {
1904  return p.first == lcRef;
1905  });
1906  if (best_cp_linked ==
1907  cPOnLayerMap[best->first].end()) // This should never happen by construction of the association maps
1908  continue;
1909  histograms.h_sharedenergy_layercl2caloparticle_vs_eta_perlayer.at(lcLayerId)->Fill(
1910  clusters[lcId].eta(), best_cp_linked->second.first / lc_en);
1911  histograms.h_sharedenergy_layercl2caloparticle_vs_phi_perlayer.at(lcLayerId)->Fill(
1912  clusters[lcId].phi(), best_cp_linked->second.first / lc_en);
1913  }
1914  } // End of loop over LayerClusters
1915 
1916  // Here Fill the plots to compute the different metrics linked to
1917  // gen-level, namely efficiency and duplicate. In this loop should restrict
1918  // only to the selected caloParaticles.
1919  for (const auto& cpId : cPSelectedIndices) {
1920  const edm::Ref<CaloParticleCollection> cpRef(caloParticleHandle, cpId);
1921  const auto& lcsIt = cPOnLayerMap.find(cpRef);
1922 
1923  std::map<unsigned int, float> cPEnergyOnLayer;
1924  for (unsigned int layerId = 0; layerId < layers * 2; ++layerId)
1925  cPEnergyOnLayer[layerId] = 0;
1926 
1927  for (const auto& simCluster : cP[cpId].simClusters()) {
1928  for (const auto& it_haf : simCluster->hits_and_fractions()) {
1929  const DetId hitid = (it_haf.first);
1930  const auto hitLayerId =
1931  recHitTools_->getLayerWithOffset(hitid) + layers * ((recHitTools_->zside(hitid) + 1) >> 1) - 1;
1932  std::unordered_map<DetId, const HGCRecHit*>::const_iterator itcheck = hitMap.find(hitid);
1933  if (itcheck != hitMap.end()) {
1934  const HGCRecHit* hit = itcheck->second;
1935  cPEnergyOnLayer[hitLayerId] += it_haf.second * hit->energy();
1936  }
1937  }
1938  }
1939 
1940  for (unsigned int layerId = 0; layerId < layers * 2; ++layerId) {
1941  if (!cPEnergyOnLayer[layerId])
1942  continue;
1943 
1944  histograms.h_denom_caloparticle_eta_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().eta());
1945  histograms.h_denom_caloparticle_phi_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().phi());
1946 
1947  if (lcsIt == cPOnLayerMap.end())
1948  continue;
1949  const auto& lcs = lcsIt->val;
1950 
1951  auto getLCLayerId = [&](const unsigned int lcId) {
1952  const auto firstHitDetId = (clusters[lcId].hitsAndFractions())[0].first;
1953  const auto lcLayerId = recHitTools_->getLayerWithOffset(firstHitDetId) +
1954  layers * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
1955  return lcLayerId;
1956  };
1957 
1958  for (const auto& lcPair : lcs) {
1959  if (getLCLayerId(lcPair.first.index()) != layerId)
1960  continue;
1961  histograms.h_score_caloparticle2layercl_perlayer.at(layerId)->Fill(lcPair.second.second);
1962  histograms.h_sharedenergy_caloparticle2layercl_perlayer.at(layerId)->Fill(
1963  lcPair.second.first / cPEnergyOnLayer[layerId], cPEnergyOnLayer[layerId]);
1964  histograms.h_energy_vs_score_caloparticle2layercl_perlayer.at(layerId)->Fill(
1965  lcPair.second.second, lcPair.second.first / cPEnergyOnLayer[layerId]);
1966  }
1967  const auto assoc = std::count_if(std::begin(lcs), std::end(lcs), [&](const auto& obj) {
1968  if (getLCLayerId(obj.first.index()) != layerId)
1969  return false;
1970  else
1971  return obj.second.second < ScoreCutCPtoLC_;
1972  });
1973  if (assoc) {
1974  histograms.h_num_caloparticle_eta_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().eta());
1975  histograms.h_num_caloparticle_phi_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().phi());
1976  if (assoc > 1) {
1977  histograms.h_numDup_caloparticle_eta_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().eta());
1978  histograms.h_numDup_caloparticle_phi_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().phi());
1979  }
1980  const auto best = std::min_element(std::begin(lcs), std::end(lcs), [&](const auto& obj1, const auto& obj2) {
1981  if (getLCLayerId(obj1.first.index()) != layerId)
1982  return false;
1983  else if (getLCLayerId(obj2.first.index()) == layerId)
1984  return obj1.second.second < obj2.second.second;
1985  else
1986  return true;
1987  });
1988  histograms.h_sharedenergy_caloparticle2layercl_vs_eta_perlayer.at(layerId)->Fill(
1989  cP[cpId].g4Tracks()[0].momentum().eta(), best->second.first / cPEnergyOnLayer[layerId]);
1990  histograms.h_sharedenergy_caloparticle2layercl_vs_phi_perlayer.at(layerId)->Fill(
1991  cP[cpId].g4Tracks()[0].momentum().phi(), best->second.first / cPEnergyOnLayer[layerId]);
1992  }
1993  }
1994  }
1995 }
1996 
1997 void HGVHistoProducerAlgo::layerClusters_to_SimClusters(
1998  const Histograms& histograms,
1999  const int count,
2000  edm::Handle<reco::CaloClusterCollection> clusterHandle,
2001  const reco::CaloClusterCollection& clusters,
2002  edm::Handle<std::vector<SimCluster>> simClusterHandle,
2003  std::vector<SimCluster> const& sC,
2004  std::vector<size_t> const& sCIndices,
2005  const std::vector<float>& mask,
2006  std::unordered_map<DetId, const HGCRecHit*> const& hitMap,
2007  unsigned int layers,
2008  const hgcal::RecoToSimCollectionWithSimClusters& scsInLayerClusterMap,
2009  const hgcal::SimToRecoCollectionWithSimClusters& lcsInSimClusterMap) const {
2010  // Here fill the plots to compute the different metrics linked to
2011  // reco-level, namely fake-rate and merge-rate. In this loop should *not*
2012  // restrict only to the selected SimClusters.
2013  for (unsigned int lcId = 0; lcId < clusters.size(); ++lcId) {
2014  if (mask[lcId] != 0.) {
2015  LogDebug("HGCalValidator") << "Skipping layer cluster " << lcId << " not belonging to mask" << std::endl;
2016  continue;
2017  }
2018  const auto firstHitDetId = (clusters[lcId].hitsAndFractions())[0].first;
2019  const auto lcLayerId =
2020  recHitTools_->getLayerWithOffset(firstHitDetId) + layers * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
2021  //Although the ones below are already created in the LC to CP association, will
2022  //recreate them here since in the post processor it looks in a specific directory.
2023  histograms.h_denom_layercl_in_simcl_eta_perlayer[count].at(lcLayerId)->Fill(clusters[lcId].eta());
2024  histograms.h_denom_layercl_in_simcl_phi_perlayer[count].at(lcLayerId)->Fill(clusters[lcId].phi());
2025  //
2026  const edm::Ref<reco::CaloClusterCollection> lcRef(clusterHandle, lcId);
2027  const auto& scsIt = scsInLayerClusterMap.find(lcRef);
2028  if (scsIt == scsInLayerClusterMap.end())
2029  continue;
2030 
2031  const auto lc_en = clusters[lcId].energy();
2032  const auto& scs = scsIt->val;
2033  // If a reconstructed LayerCluster has energy 0 but is linked to at least a
2034  // SimCluster, then his score should be 1 as set in the associator
2035  if (lc_en == 0. && !scs.empty()) {
2036  for (const auto& scPair : scs) {
2037  histograms.h_score_layercl2simcluster_perlayer[count].at(lcLayerId)->Fill(scPair.second);
2038  }
2039  continue;
2040  }
2041  //Loop through all SimClusters linked to the layer cluster under study
2042  for (const auto& scPair : scs) {
2043  LogDebug("HGCalValidator") << "layerCluster Id: \t" << lcId << "\t SC id: \t" << scPair.first.index()
2044  << "\t score \t" << scPair.second << std::endl;
2045  //This should be filled #layerClusters in layer x #linked SimClusters
2046  histograms.h_score_layercl2simcluster_perlayer[count].at(lcLayerId)->Fill(scPair.second);
2047  auto const& sc_linked =
2048  std::find_if(std::begin(lcsInSimClusterMap[scPair.first]),
2049  std::end(lcsInSimClusterMap[scPair.first]),
2050  [&lcRef](const std::pair<edm::Ref<reco::CaloClusterCollection>, std::pair<float, float>>& p) {
2051  return p.first == lcRef;
2052  });
2053  if (sc_linked ==
2054  lcsInSimClusterMap[scPair.first].end()) // This should never happen by construction of the association maps
2055  continue;
2056  histograms.h_sharedenergy_layercl2simcluster_perlayer[count].at(lcLayerId)->Fill(sc_linked->second.first / lc_en,
2057  lc_en);
2058  histograms.h_energy_vs_score_layercl2simcluster_perlayer[count].at(lcLayerId)->Fill(
2059  scPair.second, sc_linked->second.first / lc_en);
2060  }
2061  //Here he counts how many of the linked SimClusters of the layer cluster under study have a score above a certain value.
2062  const auto assoc =
2063  std::count_if(std::begin(scs), std::end(scs), [](const auto& obj) { return obj.second < ScoreCutLCtoSC_; });
2064  if (assoc) {
2065  histograms.h_num_layercl_in_simcl_eta_perlayer[count].at(lcLayerId)->Fill(clusters[lcId].eta());
2066  histograms.h_num_layercl_in_simcl_phi_perlayer[count].at(lcLayerId)->Fill(clusters[lcId].phi());
2067  if (assoc > 1) {
2068  histograms.h_numMerge_layercl_in_simcl_eta_perlayer[count].at(lcLayerId)->Fill(clusters[lcId].eta());
2069  histograms.h_numMerge_layercl_in_simcl_phi_perlayer[count].at(lcLayerId)->Fill(clusters[lcId].phi());
2070  }
2071  const auto& best = std::min_element(
2072  std::begin(scs), std::end(scs), [](const auto& obj1, const auto& obj2) { return obj1.second < obj2.second; });
2073  //From all SimClusters he founds the one with the best (lowest) score and takes his scId
2074  const auto& best_sc_linked =
2075  std::find_if(std::begin(lcsInSimClusterMap[best->first]),
2076  std::end(lcsInSimClusterMap[best->first]),
2077  [&lcRef](const std::pair<edm::Ref<reco::CaloClusterCollection>, std::pair<float, float>>& p) {
2078  return p.first == lcRef;
2079  });
2080  if (best_sc_linked ==
2081  lcsInSimClusterMap[best->first].end()) // This should never happen by construction of the association maps
2082  continue;
2083  histograms.h_sharedenergy_layercl2simcluster_vs_eta_perlayer[count].at(lcLayerId)->Fill(
2084  clusters[lcId].eta(), best_sc_linked->second.first / lc_en);
2085  histograms.h_sharedenergy_layercl2simcluster_vs_phi_perlayer[count].at(lcLayerId)->Fill(
2086  clusters[lcId].phi(), best_sc_linked->second.first / lc_en);
2087  }
2088  } // End of loop over LayerClusters
2089 
2090  // Fill the plots to compute the different metrics linked to
2091  // gen-level, namely efficiency and duplicate. In this loop should restrict
2092  // only to the selected SimClusters.
2093  for (const auto& scId : sCIndices) {
2094  const edm::Ref<SimClusterCollection> scRef(simClusterHandle, scId);
2095  const auto& lcsIt = lcsInSimClusterMap.find(scRef);
2096 
2097  std::map<unsigned int, float> sCEnergyOnLayer;
2098  for (unsigned int layerId = 0; layerId < layers * 2; ++layerId)
2099  sCEnergyOnLayer[layerId] = 0;
2100 
2101  for (const auto& it_haf : sC[scId].hits_and_fractions()) {
2102  const DetId hitid = (it_haf.first);
2103  const auto scLayerId =
2104  recHitTools_->getLayerWithOffset(hitid) + layers * ((recHitTools_->zside(hitid) + 1) >> 1) - 1;
2105  std::unordered_map<DetId, const HGCRecHit*>::const_iterator itcheck = hitMap.find(hitid);
2106  if (itcheck != hitMap.end()) {
2107  const HGCRecHit* hit = itcheck->second;
2108  sCEnergyOnLayer[scLayerId] += it_haf.second * hit->energy();
2109  }
2110  }
2111 
2112  for (unsigned int layerId = 0; layerId < layers * 2; ++layerId) {
2113  if (!sCEnergyOnLayer[layerId])
2114  continue;
2115 
2116  histograms.h_denom_simcluster_eta_perlayer[count].at(layerId)->Fill(sC[scId].eta());
2117  histograms.h_denom_simcluster_phi_perlayer[count].at(layerId)->Fill(sC[scId].phi());
2118 
2119  if (lcsIt == lcsInSimClusterMap.end())
2120  continue;
2121  const auto& lcs = lcsIt->val;
2122 
2123  auto getLCLayerId = [&](const unsigned int lcId) {
2124  const auto firstHitDetId = (clusters[lcId].hitsAndFractions())[0].first;
2125  const unsigned int lcLayerId = recHitTools_->getLayerWithOffset(firstHitDetId) +
2126  layers * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
2127  return lcLayerId;
2128  };
2129 
2130  //Loop through layer clusters linked to the SimCluster under study
2131  for (const auto& lcPair : lcs) {
2132  auto lcId = lcPair.first.index();
2133  if (mask[lcId] != 0.) {
2134  LogDebug("HGCalValidator") << "Skipping layer cluster " << lcId << " not belonging to mask" << std::endl;
2135  continue;
2136  }
2137 
2138  if (getLCLayerId(lcId) != layerId)
2139  continue;
2140  histograms.h_score_simcluster2layercl_perlayer[count].at(layerId)->Fill(lcPair.second.second);
2141  histograms.h_sharedenergy_simcluster2layercl_perlayer[count].at(layerId)->Fill(
2142  lcPair.second.first / sCEnergyOnLayer[layerId], sCEnergyOnLayer[layerId]);
2143  histograms.h_energy_vs_score_simcluster2layercl_perlayer[count].at(layerId)->Fill(
2144  lcPair.second.second, lcPair.second.first / sCEnergyOnLayer[layerId]);
2145  }
2146  const auto assoc = std::count_if(std::begin(lcs), std::end(lcs), [&](const auto& obj) {
2147  if (getLCLayerId(obj.first.index()) != layerId)
2148  return false;
2149  else
2150  return obj.second.second < ScoreCutSCtoLC_;
2151  });
2152  if (assoc) {
2153  histograms.h_num_simcluster_eta_perlayer[count].at(layerId)->Fill(sC[scId].eta());
2154  histograms.h_num_simcluster_phi_perlayer[count].at(layerId)->Fill(sC[scId].phi());
2155  if (assoc > 1) {
2156  histograms.h_numDup_simcluster_eta_perlayer[count].at(layerId)->Fill(sC[scId].eta());
2157  histograms.h_numDup_simcluster_phi_perlayer[count].at(layerId)->Fill(sC[scId].phi());
2158  }
2159  const auto best = std::min_element(std::begin(lcs), std::end(lcs), [&](const auto& obj1, const auto& obj2) {
2160  if (getLCLayerId(obj1.first.index()) != layerId)
2161  return false;
2162  else if (getLCLayerId(obj2.first.index()) == layerId)
2163  return obj1.second.second < obj2.second.second;
2164  else
2165  return true;
2166  });
2167  histograms.h_sharedenergy_simcluster2layercl_vs_eta_perlayer[count].at(layerId)->Fill(
2168  sC[scId].eta(), best->second.first / sCEnergyOnLayer[layerId]);
2169  histograms.h_sharedenergy_simcluster2layercl_vs_phi_perlayer[count].at(layerId)->Fill(
2170  sC[scId].phi(), best->second.first / sCEnergyOnLayer[layerId]);
2171  }
2172  }
2173  }
2174 }
2175 
2176 void HGVHistoProducerAlgo::fill_generic_cluster_histos(const Histograms& histograms,
2177  const int count,
2178  edm::Handle<reco::CaloClusterCollection> clusterHandle,
2179  const reco::CaloClusterCollection& clusters,
2180  const Density& densities,
2181  edm::Handle<std::vector<CaloParticle>> caloParticleHandle,
2182  std::vector<CaloParticle> const& cP,
2183  std::vector<size_t> const& cPIndices,
2184  std::vector<size_t> const& cPSelectedIndices,
2185  std::unordered_map<DetId, const HGCRecHit*> const& hitMap,
2186  std::map<double, double> cummatbudg,
2187  unsigned int layers,
2188  std::vector<int> thicknesses,
2189  const hgcal::RecoToSimCollection& cpsInLayerClusterMap,
2190  const hgcal::SimToRecoCollection& cPOnLayerMap) const {
2191  //Each event to be treated as two events: an event in +ve endcap,
2192  //plus another event in -ve endcap. In this spirit there will be
2193  //a layer variable (layerid) that maps the layers in :
2194  //-z: 0->51
2195  //+z: 52->103
2196 
2197  //To keep track of total num of layer clusters per layer
2198  //tnlcpl[layerid]
2199  std::vector<int> tnlcpl(1000, 0); //tnlcpl.clear(); tnlcpl.reserve(1000);
2200 
2201  //To keep track of the total num of clusters per thickness in plus and in minus endcaps
2202  std::map<std::string, int> tnlcpthplus;
2203  tnlcpthplus.clear();
2204  std::map<std::string, int> tnlcpthminus;
2205  tnlcpthminus.clear();
2206  //At the beginning of the event all layers should be initialized to zero total clusters per thickness
2207  for (std::vector<int>::iterator it = thicknesses.begin(); it != thicknesses.end(); ++it) {
2208  tnlcpthplus.insert(std::pair<std::string, int>(std::to_string(*it), 0));
2209  tnlcpthminus.insert(std::pair<std::string, int>(std::to_string(*it), 0));
2210  }
2211  //To keep track of the total num of clusters with mixed thickness hits per event
2212  tnlcpthplus.insert(std::pair<std::string, int>("mixed", 0));
2213  tnlcpthminus.insert(std::pair<std::string, int>("mixed", 0));
2214 
2215  layerClusters_to_CaloParticles(histograms,
2216  clusterHandle,
2217  clusters,
2218  caloParticleHandle,
2219  cP,
2220  cPIndices,
2221  cPSelectedIndices,
2222  hitMap,
2223  layers,
2224  cpsInLayerClusterMap,
2225  cPOnLayerMap);
2226 
2227  //To find out the total amount of energy clustered per layer
2228  //Initialize with zeros because I see clear gives weird numbers.
2229  std::vector<double> tecpl(1000, 0.0); //tecpl.clear(); tecpl.reserve(1000);
2230  //for the longitudinal depth barycenter
2231  std::vector<double> ldbar(1000, 0.0); //ldbar.clear(); ldbar.reserve(1000);
2232 
2233  // Need to compare with the total amount of energy coming from CaloParticles
2234  double caloparteneplus = 0.;
2235  double caloparteneminus = 0.;
2236  for (const auto& cpId : cPIndices) {
2237  if (cP[cpId].eta() >= 0.) {
2238  caloparteneplus = caloparteneplus + cP[cpId].energy();
2239  } else if (cP[cpId].eta() < 0.) {
2240  caloparteneminus = caloparteneminus + cP[cpId].energy();
2241  }
2242  }
2243 
2244  // loop through clusters of the event
2245  for (const auto& lcId : clusters) {
2246  const auto seedid = lcId.seed();
2247  const double seedx = recHitTools_->getPosition(seedid).x();
2248  const double seedy = recHitTools_->getPosition(seedid).y();
2249  DetId maxid = findmaxhit(lcId, hitMap);
2250 
2251  // const DetId maxid = lcId.max();
2252  double maxx = recHitTools_->getPosition(maxid).x();
2253  double maxy = recHitTools_->getPosition(maxid).y();
2254 
2255  //Auxillary variables to count the number of different kind of hits in each cluster
2256  int nthhits120p = 0;
2257  int nthhits200p = 0;
2258  int nthhits300p = 0;
2259  int nthhitsscintp = 0;
2260  int nthhits120m = 0;
2261  int nthhits200m = 0;
2262  int nthhits300m = 0;
2263  int nthhitsscintm = 0;
2264  //For the hits thickness of the layer cluster.
2265  double thickness = 0.;
2266  //The layer the cluster belongs to. As mentioned in the mapping above, it takes into account -z and +z.
2267  int layerid = 0;
2268  // Need another layer variable for the longitudinal material budget file reading.
2269  //In this case need no distinction between -z and +z.
2270  int lay = 0;
2271  // Need to save the combination thick_lay
2272  std::string istr = "";
2273  //boolean to check for the layer that the cluster belong to. Maybe later will check all the layer hits.
2274  bool cluslay = true;
2275  //zside that the current cluster belongs to.
2276  int zside = 0;
2277 
2278  const auto& hits_and_fractions = lcId.hitsAndFractions();
2279  for (std::vector<std::pair<DetId, float>>::const_iterator it_haf = hits_and_fractions.begin();
2280  it_haf != hits_and_fractions.end();
2281  ++it_haf) {
2282  const DetId rh_detid = it_haf->first;
2283  //The layer that the current hit belongs to
2284  layerid = recHitTools_->getLayerWithOffset(rh_detid) + layers * ((recHitTools_->zside(rh_detid) + 1) >> 1) - 1;
2285  lay = recHitTools_->getLayerWithOffset(rh_detid);
2286  zside = recHitTools_->zside(rh_detid);
2287  if (rh_detid.det() == DetId::Forward || rh_detid.det() == DetId::HGCalEE || rh_detid.det() == DetId::HGCalHSi) {
2288  thickness = recHitTools_->getSiThickness(rh_detid);
2289  } else if (rh_detid.det() == DetId::HGCalHSc) {
2290  thickness = -1;
2291  } else {
2292  LogDebug("HGCalValidator") << "These are HGCal layer clusters, you shouldn't be here !!! " << layerid << "\n";
2293  continue;
2294  }
2295 
2296  //Count here only once the layer cluster and save the combination thick_layerid
2297  std::string curistr = std::to_string((int)thickness);
2298  std::string lay_string = std::to_string(layerid);
2299  while (lay_string.size() < 2)
2300  lay_string.insert(0, "0");
2301  curistr += "_" + lay_string;
2302  if (cluslay) {
2303  tnlcpl[layerid]++;
2304  istr = curistr;
2305  cluslay = false;
2306  }
2307 
2308  if ((thickness == 120.) && (recHitTools_->zside(rh_detid) > 0.)) {
2309  nthhits120p++;
2310  } else if ((thickness == 120.) && (recHitTools_->zside(rh_detid) < 0.)) {
2311  nthhits120m++;
2312  } else if ((thickness == 200.) && (recHitTools_->zside(rh_detid) > 0.)) {
2313  nthhits200p++;
2314  } else if ((thickness == 200.) && (recHitTools_->zside(rh_detid) < 0.)) {
2315  nthhits200m++;
2316  } else if ((thickness == 300.) && (recHitTools_->zside(rh_detid) > 0.)) {
2317  nthhits300p++;
2318  } else if ((thickness == 300.) && (recHitTools_->zside(rh_detid) < 0.)) {
2319  nthhits300m++;
2320  } else if ((thickness == -1) && (recHitTools_->zside(rh_detid) > 0.)) {
2321  nthhitsscintp++;
2322  } else if ((thickness == -1) && (recHitTools_->zside(rh_detid) < 0.)) {
2323  nthhitsscintm++;
2324  } else { //assert(0);
2325  LogDebug("HGCalValidator")
2326  << " You are running a geometry that contains thicknesses different than the normal ones. "
2327  << "\n";
2328  }
2329 
2330  std::unordered_map<DetId, const HGCRecHit*>::const_iterator itcheck = hitMap.find(rh_detid);
2331  if (itcheck == hitMap.end()) {
2332  std::ostringstream st1;
2333  if ((rh_detid.det() == DetId::HGCalEE) || (rh_detid.det() == DetId::HGCalHSi)) {
2334  st1 << HGCSiliconDetId(rh_detid);
2335  } else if (rh_detid.det() == DetId::HGCalHSc) {
2336  st1 << HGCScintillatorDetId(rh_detid);
2337  } else {
2338  st1 << HFNoseDetId(rh_detid);
2339  }
2340  LogDebug("HGCalValidator") << " You shouldn't be here - Unable to find a hit " << rh_detid.rawId() << " "
2341  << rh_detid.det() << " " << st1.str() << "\n";
2342  continue;
2343  }
2344 
2345  const HGCRecHit* hit = itcheck->second;
2346 
2347  //Here for the per cell plots
2348  //----
2349  const double hit_x = recHitTools_->getPosition(rh_detid).x();
2350  const double hit_y = recHitTools_->getPosition(rh_detid).y();
2351  double distancetoseed = distance(seedx, seedy, hit_x, hit_y);
2352  double distancetomax = distance(maxx, maxy, hit_x, hit_y);
2353  if (distancetoseed != 0. && histograms.h_distancetoseedcell_perthickperlayer.count(curistr)) {
2354  histograms.h_distancetoseedcell_perthickperlayer.at(curistr)->Fill(distancetoseed);
2355  }
2356  //----
2357  if (distancetoseed != 0. && histograms.h_distancetoseedcell_perthickperlayer_eneweighted.count(curistr)) {
2358  histograms.h_distancetoseedcell_perthickperlayer_eneweighted.at(curistr)->Fill(distancetoseed, hit->energy());
2359  }
2360  //----
2361  if (distancetomax != 0. && histograms.h_distancetomaxcell_perthickperlayer.count(curistr)) {
2362  histograms.h_distancetomaxcell_perthickperlayer.at(curistr)->Fill(distancetomax);
2363  }
2364  //----
2365  if (distancetomax != 0. && histograms.h_distancetomaxcell_perthickperlayer_eneweighted.count(curistr)) {
2366  histograms.h_distancetomaxcell_perthickperlayer_eneweighted.at(curistr)->Fill(distancetomax, hit->energy());
2367  }
2368 
2369  //Let's check the density
2370  std::map<DetId, float>::const_iterator dit = densities.find(rh_detid);
2371  if (dit == densities.end()) {
2372  std::ostringstream st1;
2373  if ((rh_detid.det() == DetId::HGCalEE) || (rh_detid.det() == DetId::HGCalHSi)) {
2374  st1 << HGCSiliconDetId(rh_detid);
2375  } else if (rh_detid.det() == DetId::HGCalHSc) {
2376  st1 << HGCScintillatorDetId(rh_detid);
2377  } else {
2378  st1 << HFNoseDetId(rh_detid);
2379  }
2380  LogDebug("HGCalValidator") << " You shouldn't be here - Unable to find a density " << rh_detid.rawId() << " "
2381  << rh_detid.det() << " " << st1.str() << "\n";
2382  continue;
2383  }
2384 
2385  if (histograms.h_cellsenedens_perthick.count((int)thickness)) {
2386  histograms.h_cellsenedens_perthick.at((int)thickness)->Fill(dit->second);
2387  }
2388 
2389  } // end of loop through hits and fractions
2390 
2391  //Check for simultaneously having hits of different kind. Checking at least two combinations is sufficient.
2392  if ((nthhits120p != 0 && nthhits200p != 0) || (nthhits120p != 0 && nthhits300p != 0) ||
2393  (nthhits120p != 0 && nthhitsscintp != 0) || (nthhits200p != 0 && nthhits300p != 0) ||
2394  (nthhits200p != 0 && nthhitsscintp != 0) || (nthhits300p != 0 && nthhitsscintp != 0)) {
2395  tnlcpthplus["mixed"]++;
2396  } else if ((nthhits120p != 0 || nthhits200p != 0 || nthhits300p != 0 || nthhitsscintp != 0)) {
2397  //This is a cluster with hits of one kind
2398  tnlcpthplus[std::to_string((int)thickness)]++;
2399  }
2400  if ((nthhits120m != 0 && nthhits200m != 0) || (nthhits120m != 0 && nthhits300m != 0) ||
2401  (nthhits120m != 0 && nthhitsscintm != 0) || (nthhits200m != 0 && nthhits300m != 0) ||
2402  (nthhits200m != 0 && nthhitsscintm != 0) || (nthhits300m != 0 && nthhitsscintm != 0)) {
2403  tnlcpthminus["mixed"]++;
2404  } else if ((nthhits120m != 0 || nthhits200m != 0 || nthhits300m != 0 || nthhitsscintm != 0)) {
2405  //This is a cluster with hits of one kind
2406  tnlcpthminus[std::to_string((int)thickness)]++;
2407  }
2408 
2409  //To find the thickness with the biggest amount of cells
2410  std::vector<int> bigamoth;
2411  bigamoth.clear();
2412  if (zside > 0) {
2413  bigamoth.push_back(nthhits120p);
2414  bigamoth.push_back(nthhits200p);
2415  bigamoth.push_back(nthhits300p);
2416  bigamoth.push_back(nthhitsscintp);
2417  } else if (zside < 0) {
2418  bigamoth.push_back(nthhits120m);
2419  bigamoth.push_back(nthhits200m);
2420  bigamoth.push_back(nthhits300m);
2421  bigamoth.push_back(nthhitsscintm);
2422  }
2423  auto bgth = std::max_element(bigamoth.begin(), bigamoth.end());
2424  istr = std::to_string(thicknesses[std::distance(bigamoth.begin(), bgth)]);
2425  std::string lay_string = std::to_string(layerid);
2426  while (lay_string.size() < 2)
2427  lay_string.insert(0, "0");
2428  istr += "_" + lay_string;
2429 
2430  //Here for the per cluster plots that need the thickness_layer info
2431  if (histograms.h_cellsnum_perthickperlayer.count(istr)) {
2432  histograms.h_cellsnum_perthickperlayer.at(istr)->Fill(hits_and_fractions.size());
2433  }
2434 
2435  //Now, with the distance between seed and max cell.
2436  double distancebetseedandmax = distance(seedx, seedy, maxx, maxy);
2437  //The thickness_layer combination in this case will use the thickness of the seed as a convention.
2438  std::string seedstr = std::to_string((int)recHitTools_->getSiThickness(seedid)) + "_" + std::to_string(layerid);
2439  seedstr += "_" + lay_string;
2440  if (histograms.h_distancebetseedandmaxcell_perthickperlayer.count(seedstr)) {
2441  histograms.h_distancebetseedandmaxcell_perthickperlayer.at(seedstr)->Fill(distancebetseedandmax);
2442  }
2443  const auto lc_en = lcId.energy();
2444  if (histograms.h_distancebetseedandmaxcellvsclusterenergy_perthickperlayer.count(seedstr)) {
2445  histograms.h_distancebetseedandmaxcellvsclusterenergy_perthickperlayer.at(seedstr)->Fill(distancebetseedandmax,
2446  lc_en);
2447  }
2448 
2449  //Energy clustered per layer
2450  tecpl[layerid] = tecpl[layerid] + lc_en;
2451  ldbar[layerid] = ldbar[layerid] + lc_en * cummatbudg[(double)lay];
2452 
2453  } //end of loop through clusters of the event
2454 
2455  // First a couple of variables to keep the sum of the energy of all clusters
2456  double sumeneallcluspl = 0.;
2457  double sumeneallclusmi = 0.;
2458  // and the longitudinal variable
2459  double sumldbarpl = 0.;
2460  double sumldbarmi = 0.;
2461  //Per layer : Loop 0->103
2462  for (unsigned ilayer = 0; ilayer < layers * 2; ++ilayer) {
2463  if (histograms.h_clusternum_perlayer.count(ilayer)) {
2464  histograms.h_clusternum_perlayer.at(ilayer)->Fill(tnlcpl[ilayer]);
2465  }
2466  // Two times one for plus and one for minus
2467  //First with the -z endcap
2468  if (ilayer < layers) {
2469  if (histograms.h_energyclustered_perlayer.count(ilayer)) {
2470  if (caloparteneminus != 0.) {
2471  histograms.h_energyclustered_perlayer.at(ilayer)->Fill(100. * tecpl[ilayer] / caloparteneminus);
2472  }
2473  }
2474  //Keep here the total energy for the event in -z
2475  sumeneallclusmi = sumeneallclusmi + tecpl[ilayer];
2476  //And for the longitudinal variable
2477  sumldbarmi = sumldbarmi + ldbar[ilayer];
2478  } else { //Then for the +z
2479  if (histograms.h_energyclustered_perlayer.count(ilayer)) {
2480  if (caloparteneplus != 0.) {
2481  histograms.h_energyclustered_perlayer.at(ilayer)->Fill(100. * tecpl[ilayer] / caloparteneplus);
2482  }
2483  }
2484  //Keep here the total energy for the event in -z
2485  sumeneallcluspl = sumeneallcluspl + tecpl[ilayer];
2486  //And for the longitudinal variable
2487  sumldbarpl = sumldbarpl + ldbar[ilayer];
2488  } //end of +z loop
2489 
2490  } //end of loop over layers
2491 
2492  //Per thickness
2493  for (std::vector<int>::iterator it = thicknesses.begin(); it != thicknesses.end(); ++it) {
2494  if (histograms.h_clusternum_perthick.count(*it)) {
2495  histograms.h_clusternum_perthick.at(*it)->Fill(tnlcpthplus[std::to_string(*it)]);
2496  histograms.h_clusternum_perthick.at(*it)->Fill(tnlcpthminus[std::to_string(*it)]);
2497  }
2498  }
2499  //Mixed thickness clusters
2500  histograms.h_mixedhitscluster_zplus[count]->Fill(tnlcpthplus["mixed"]);
2501  histograms.h_mixedhitscluster_zminus[count]->Fill(tnlcpthminus["mixed"]);
2502 
2503  //Total energy clustered from all layer clusters (fraction)
2504  if (caloparteneplus != 0.) {
2505  histograms.h_energyclustered_zplus[count]->Fill(100. * sumeneallcluspl / caloparteneplus);
2506  }
2507  if (caloparteneminus != 0.) {
2508  histograms.h_energyclustered_zminus[count]->Fill(100. * sumeneallclusmi / caloparteneminus);
2509  }
2510 
2511  //For the longitudinal depth barycenter
2512  histograms.h_longdepthbarycentre_zplus[count]->Fill(sumldbarpl / sumeneallcluspl);
2513  histograms.h_longdepthbarycentre_zminus[count]->Fill(sumldbarmi / sumeneallclusmi);
2514 }
2515 
2516 void HGVHistoProducerAlgo::tracksters_to_SimTracksters(
2517  const Histograms& histograms,
2518  const int count,
2519  const ticl::TracksterCollection& tracksters,
2520  const reco::CaloClusterCollection& layerClusters,
2521  const ticl::TracksterCollection& simTSs,
2522  const validationType valType,
2523  const ticl::TracksterCollection& simTSs_fromCP,
2524  const std::map<unsigned int, std::vector<unsigned int>>& cpToSc_SimTrackstersMap,
2525  std::vector<SimCluster> const& sC,
2526  const edm::ProductID& cPHandle_id,
2527  std::vector<CaloParticle> const& cP,
2528  std::vector<size_t> const& cPIndices,
2529  std::vector<size_t> const& cPSelectedIndices,
2530  std::unordered_map<DetId, const HGCRecHit*> const& hitMap,
2531  unsigned int layers) const {
2532  const auto nTracksters = tracksters.size();
2533  const auto nSimTracksters = simTSs.size();
2534 
2535  std::unordered_map<DetId, std::vector<HGVHistoProducerAlgo::detIdInfoInCluster>> detIdSimTSId_Map;
2536  std::unordered_map<DetId, std::vector<HGVHistoProducerAlgo::detIdInfoInTrackster>> detIdToTracksterId_Map;
2537  std::vector<int> tracksters_FakeMerge(nTracksters, 0);
2538  std::vector<int> tracksters_PurityDuplicate(nTracksters, 0);
2539 
2540  // This vector contains the ids of the SimTracksters contributing with at least one hit to the Trackster and the reconstruction error
2541  //stsInTrackster[trackster][STSids]
2542  //Connects a Trackster with all related SimTracksters.
2543  std::vector<std::vector<std::pair<unsigned int, float>>> stsInTrackster;
2544  stsInTrackster.resize(nTracksters);
2545 
2546  // cPOnLayer[caloparticle]:
2547  //1. the sum of all rechits energy times fraction of the relevant simhit related to that calo particle.
2548  //2. the hits and fractions of that calo particle.
2549  //3. the layer clusters with matched rechit id.
2550  std::unordered_map<int, caloParticleOnLayer> cPOnLayer;
2551  std::unordered_map<int, std::vector<caloParticleOnLayer>> sCOnLayer;
2552  //Consider CaloParticles coming from the hard scatterer, excluding the PU contribution.
2553  for (const auto cpIndex : cPIndices) {
2554  cPOnLayer[cpIndex].caloParticleId = cpIndex;
2555  cPOnLayer[cpIndex].energy = 0.f;
2556  cPOnLayer[cpIndex].hits_and_fractions.clear();
2557  const auto nSC_inCP = sC.size();
2558  sCOnLayer[cpIndex].resize(nSC_inCP);
2559  for (unsigned int iSC = 0; iSC < nSC_inCP; iSC++) {
2560  sCOnLayer[cpIndex][iSC].caloParticleId = cpIndex;
2561  sCOnLayer[cpIndex][iSC].energy = 0.f;
2562  sCOnLayer[cpIndex][iSC].hits_and_fractions.clear();
2563  }
2564  }
2565 
2566  auto getCPId = [](const ticl::Trackster& simTS,
2567  const unsigned int iSTS,
2568  const edm::ProductID& cPHandle_id,
2569  const std::map<unsigned int, std::vector<unsigned int>>& cpToSc_SimTrackstersMap,
2570  const ticl::TracksterCollection& simTSs_fromCP) {
2571  unsigned int cpId = -1;
2572 
2573  const auto productID = simTS.seedID();
2574  if (productID == cPHandle_id) {
2575  cpId = simTS.seedIndex();
2576  } else { // SimTrackster from SimCluster
2577  const auto findSimTSFromCP = std::find_if(
2578  std::begin(cpToSc_SimTrackstersMap),
2579  std::end(cpToSc_SimTrackstersMap),
2580  [&](const std::pair<unsigned int, std::vector<unsigned int>>& cpToScs) {
2581  return std::find(std::begin(cpToScs.second), std::end(cpToScs.second), iSTS) != std::end(cpToScs.second);
2582  });
2583  if (findSimTSFromCP != std::end(cpToSc_SimTrackstersMap)) {
2584  cpId = simTSs_fromCP[findSimTSFromCP->first].seedIndex();
2585  }
2586  }
2587 
2588  return cpId;
2589  };
2590 
2591  auto getLCId = [](const std::vector<unsigned int>& tst_vertices,
2592  const reco::CaloClusterCollection& layerClusters,
2593  const DetId& hitid) {
2594  unsigned int lcId = -1;
2595  std::for_each(std::begin(tst_vertices), std::end(tst_vertices), [&](unsigned int idx) {
2596  const auto& lc_haf = layerClusters[idx].hitsAndFractions();
2597  const auto& hitFound = std::find_if(std::begin(lc_haf),
2598  std::end(lc_haf),
2599  [&hitid](const std::pair<DetId, float>& v) { return v.first == hitid; });
2600  if (hitFound != lc_haf.end()) // not all hits may be clusterized
2601  lcId = idx;
2602  });
2603  return lcId;
2604  };
2605 
2606  for (unsigned int iSTS = 0; iSTS < nSimTracksters; ++iSTS) {
2607  const auto cpId = getCPId(simTSs[iSTS], iSTS, cPHandle_id, cpToSc_SimTrackstersMap, simTSs_fromCP);
2608  if (std::find(cPIndices.begin(), cPIndices.end(), cpId) == cPIndices.end())
2609  continue;
2610 
2611  // Loop through SimClusters
2612  for (const auto& simCluster : cP[cpId].simClusters()) {
2613  auto iSim = simTSs[iSTS].seedIndex();
2614  if (simTSs[iSTS].seedID() != cPHandle_id) { // SimTrackster from SimCluster
2615  if (iSim != (&(*simCluster) - &(sC[0])))
2616  continue;
2617  } else
2618  iSim = 0;
2619 
2620  for (const auto& it_haf : simCluster->hits_and_fractions()) {
2621  const auto hitid = (it_haf.first);
2622  const auto lcId = getLCId(simTSs[iSTS].vertices(), layerClusters, hitid);
2623  //V9:maps the layers in -z: 0->51 and in +z: 52->103
2624  //V10:maps the layers in -z: 0->49 and in +z: 50->99
2625  const auto itcheck = hitMap.find(hitid);
2626  //Check whether the current hit belonging to sim cluster has a reconstructed hit
2627  if ((valType == 0 && itcheck != hitMap.end()) || (valType > 0 && int(lcId) >= 0)) {
2628  float lcFraction = 0;
2629  if (valType > 0) {
2630  const auto iLC = std::find(simTSs[iSTS].vertices().begin(), simTSs[iSTS].vertices().end(), lcId);
2631  lcFraction =
2632  1.f / simTSs[iSTS].vertex_multiplicity(std::distance(std::begin(simTSs[iSTS].vertices()), iLC));
2633  }
2634  const auto elemId = (valType == 0) ? hitid : lcId;
2635  const auto elemFr = (valType == 0) ? it_haf.second : lcFraction;
2636  //Since the current hit from SimCluster has a reconstructed hit {with the same detid, belonging to the corresponding SimTrackster},
2637  //make a map that will connect a {detid,lcID} with:
2638  //1. the SimTracksters that have a SimCluster with {sim hits in, LCs containing} that cell via SimTrackster id.
2639  //2. the sum of all {SimHits, LCs} fractions that contributes to that detid.
2640  //So, keep in mind that in case of multiple CaloParticles contributing in the same cell
2641  //the fraction is the sum over all calo particles. So, something like:
2642  //detid: (caloparticle 1, sum of hits fractions in that detid over all cp) , (caloparticle 2, sum of hits fractions in that detid over all cp), (caloparticle 3, sum of hits fractions in that detid over all cp) ...
2643  if (detIdSimTSId_Map.find(elemId) == detIdSimTSId_Map.end()) {
2644  detIdSimTSId_Map[elemId] = std::vector<HGVHistoProducerAlgo::detIdInfoInCluster>();
2645  detIdSimTSId_Map[elemId].emplace_back(HGVHistoProducerAlgo::detIdInfoInCluster{iSTS, elemFr});
2646  } else {
2647  auto findSTSIt =
2648  std::find(detIdSimTSId_Map[elemId].begin(),
2649  detIdSimTSId_Map[elemId].end(),
2650  HGVHistoProducerAlgo::detIdInfoInCluster{
2651  iSTS, 0}); // only the first element is used for the matching (overloaded operator==)
2652  if (findSTSIt != detIdSimTSId_Map[elemId].end()) {
2653  if (valType == 0)
2654  findSTSIt->fraction += elemFr;
2655  } else {
2656  detIdSimTSId_Map[elemId].emplace_back(HGVHistoProducerAlgo::detIdInfoInCluster{iSTS, elemFr});
2657  }
2658  }
2659  const auto hitEn = itcheck->second->energy();
2660  //Since the current hit from sim cluster has a reconstructed hit with the same detid,
2661  //fill the cPOnLayer[caloparticle][layer] object with energy (sum of all rechits energy times fraction
2662  //of the relevant simhit) and keep the hit (detid and fraction) that contributed.
2663  cPOnLayer[cpId].energy += it_haf.second * hitEn;
2664  sCOnLayer[cpId][iSim].energy += elemFr * hitEn;
2665  // Need to compress the hits and fractions in order to have a
2666  // reasonable score between CP and LC. Imagine, for example, that a
2667  // CP has detID X used by 2 SimClusters with different fractions. If
2668  // a single LC uses X with fraction 1 and is compared to the 2
2669  // contributions separately, it will be assigned a score != 0, which
2670  // is wrong.
2671  auto& haf = cPOnLayer[cpId].hits_and_fractions;
2672  auto found = std::find_if(
2673  std::begin(haf), std::end(haf), [&hitid](const std::pair<DetId, float>& v) { return v.first == hitid; });
2674  if (found != haf.end())
2675  found->second += it_haf.second;
2676  else
2677  haf.emplace_back(hitid, it_haf.second);
2678  // Same for sCOnLayer
2679  auto& haf_sc = sCOnLayer[cpId][iSim].hits_and_fractions;
2680  auto found_sc = std::find_if(std::begin(haf_sc),
2681  std::end(haf_sc),
2682  [&hitid](const std::pair<DetId, float>& v) { return v.first == hitid; });
2683  if (found_sc != haf_sc.end())
2684  found_sc->second += it_haf.second;
2685  else
2686  haf_sc.emplace_back(hitid, it_haf.second);
2687  }
2688  } // end of loop through SimHits
2689  } // end of loop through SimClusters
2690  } // end of loop through SimTracksters
2691 
2692  auto apply_LCMultiplicity = [](const ticl::Trackster& trackster, const reco::CaloClusterCollection& layerClusters) {
2693  std::vector<std::pair<DetId, float>> hits_and_fractions_norm;
2694  int lcInTst = 0;
2695  std::for_each(std::begin(trackster.vertices()), std::end(trackster.vertices()), [&](unsigned int idx) {
2696  const auto fraction = 1.f / trackster.vertex_multiplicity(lcInTst++);
2697  for (const auto& cell : layerClusters[idx].hitsAndFractions()) {
2698  hits_and_fractions_norm.emplace_back(
2699  cell.first, cell.second * fraction); // cell.second is the hit fraction in the layerCluster
2700  }
2701  });
2702  return hits_and_fractions_norm;
2703  };
2704 
2705  auto ScoreCutSTStoTSPurDup = ScoreCutSTStoTSPurDup_[0];
2706  auto ScoreCutTStoSTSFakeMerge = ScoreCutTStoSTSFakeMerge_[0];
2707  // Loop through Tracksters
2708  for (unsigned int tstId = 0; tstId < nTracksters; ++tstId) {
2709  const auto& tst = tracksters[tstId];
2710  if (tstId == 0)
2711  if ((valType > 0) && (tst.ticlIteration() == ticl::Trackster::SIM)) {
2712  ScoreCutSTStoTSPurDup = ScoreCutSTStoTSPurDup_[valType];
2713  ScoreCutTStoSTSFakeMerge = ScoreCutTStoSTSFakeMerge_[valType];
2714  }
2715 
2716  if (tst.vertices().empty())
2717  continue;
2718 
2719  std::unordered_map<unsigned, float> CPEnergyInTS;
2720  int maxCPId_byNumberOfHits = -1;
2721  unsigned int maxCPNumberOfHitsInTS = 0;
2722  int maxCPId_byEnergy = -1;
2723  float maxEnergySharedTSandCP = 0.f;
2724  float energyFractionOfTSinCP = 0.f;
2725  float energyFractionOfCPinTS = 0.f;
2726 
2727  //In case of matched rechit-simhit, so matched
2728  //CaloParticle-LayerCluster-Trackster, he counts and saves the number of
2729  //rechits related to the maximum energy CaloParticle out of all
2730  //CaloParticles related to that layer cluster and Trackster.
2731 
2732  std::unordered_map<unsigned, unsigned> occurrencesCPinTS;
2733  unsigned int numberOfNoiseHitsInTS = 0;
2734  unsigned int numberOfHaloHitsInTS = 0;
2735 
2736  const auto tst_hitsAndFractions = apply_LCMultiplicity(tst, layerClusters);
2737  const auto numberOfHitsInTS = tst_hitsAndFractions.size();
2738 
2739  //hitsToCaloParticleId is a vector of ints, one for each rechit of the
2740  //layer cluster under study. If negative, there is no simhit from any CaloParticle related.
2741  //If positive, at least one CaloParticle has been found with matched simhit.
2742  //In more detail:
2743  // 1. hitsToCaloParticleId[iHit] = -3
2744  // TN: These represent Halo Cells(N) that have not been
2745  // assigned to any CaloParticle (hence the T).
2746  // 2. hitsToCaloParticleId[iHit] = -2
2747  // FN: There represent Halo Cells(N) that have been assigned
2748  // to a CaloParticle (hence the F, since those should have not been marked as halo)
2749  // 3. hitsToCaloParticleId[iHit] = -1
2750  // FP: These represent Real Cells(P) that have not been
2751  // assigned to any CaloParticle (hence the F, since these are fakes)
2752  // 4. hitsToCaloParticleId[iHit] >= 0
2753  // TP There represent Real Cells(P) that have been assigned
2754  // to a CaloParticle (hence the T)
2755  std::vector<int> hitsToCaloParticleId(numberOfHitsInTS);
2756 
2757  //Loop through the hits of the trackster under study
2758  for (unsigned int iHit = 0; iHit < numberOfHitsInTS; iHit++) {
2759  const auto rh_detid = tst_hitsAndFractions[iHit].first;
2760  const auto rhFraction = tst_hitsAndFractions[iHit].second;
2761 
2762  const auto lcId_r = getLCId(tst.vertices(), layerClusters, rh_detid);
2763  const auto iLC_r = std::find(tst.vertices().begin(), tst.vertices().end(), lcId_r);
2764  const auto lcFraction_r = 1.f / tst.vertex_multiplicity(std::distance(std::begin(tst.vertices()), iLC_r));
2765 
2766  //Make a map that will connect a detid (that belongs to a rechit of the layer cluster under study,
2767  //no need to save others) with:
2768  //1. the layer clusters that have rechits in that detid
2769  //2. the fraction of the rechit of each layer cluster that contributes to that detid.
2770  //So, something like:
2771  //detid: (layer cluster 1, hit fraction) , (layer cluster 2, hit fraction), (layer cluster 3, hit fraction) ...
2772  //here comparing with the calo particle map above the
2773  if (detIdToTracksterId_Map.find(rh_detid) == detIdToTracksterId_Map.end()) {
2774  detIdToTracksterId_Map[rh_detid] = std::vector<HGVHistoProducerAlgo::detIdInfoInTrackster>();
2775  detIdToTracksterId_Map[rh_detid].emplace_back(
2776  HGVHistoProducerAlgo::detIdInfoInTrackster{tstId, lcId_r, rhFraction});
2777  } else {
2778  auto findTSIt =
2779  std::find(detIdToTracksterId_Map[rh_detid].begin(),
2780  detIdToTracksterId_Map[rh_detid].end(),
2781  HGVHistoProducerAlgo::detIdInfoInTrackster{
2782  tstId, 0, 0}); // only the first element is used for the matching (overloaded operator==)
2783  if (findTSIt != detIdToTracksterId_Map[rh_detid].end()) {
2784  if (valType == 0)
2785  findTSIt->fraction += rhFraction;
2786  } else {
2787  detIdToTracksterId_Map[rh_detid].emplace_back(
2788  HGVHistoProducerAlgo::detIdInfoInTrackster{tstId, lcId_r, rhFraction});
2789  }
2790  }
2791 
2792  // if the fraction is zero or the hit does not belong to any calo
2793  // particle, set the caloparticleId for the hit to -1 this will
2794  // contribute to the number of noise hits
2795  // MR Remove the case in which the fraction is 0, since this could be a
2796  // real hit that has been marked as halo.
2797  if (rhFraction == 0.) {
2798  hitsToCaloParticleId[iHit] = -2;
2799  numberOfHaloHitsInTS++;
2800  }
2801 
2802  // Check whether the RecHit of the trackster under study has a SimHit in the same cell
2803  const auto elemId = (valType == 0) ? rh_detid.rawId() : lcId_r;
2804  const auto recoFr = (valType == 0) ? rhFraction : lcFraction_r;
2805  const auto& hit_find_in_STS = detIdSimTSId_Map.find(elemId);
2806  if (hit_find_in_STS == detIdSimTSId_Map.end()) {
2807  hitsToCaloParticleId[iHit] -= 1;
2808  } else {
2809  // Since the hit is belonging to the layer cluster, it must be also in the rechits map
2810  const auto hitEn = hitMap.find(rh_detid)->second->energy();
2811  //const auto layerId =
2812  //recHitTools_->getLayerWithOffset(rh_detid) + layers * ((recHitTools_->zside(rh_detid) + 1) >> 1) - 1;
2813  //0;
2814 
2815  auto maxCPEnergyInTS = 0.f;
2816  auto maxCPId = -1;
2817  for (const auto& h : hit_find_in_STS->second) {
2818  const auto shared_fraction = std::min(recoFr, h.fraction);
2819  const auto iSTS = h.clusterId;
2820  const auto& simTS = simTSs[iSTS];
2821  auto iSim = simTS.seedIndex();
2822  if (simTSs[iSTS].seedID() == cPHandle_id) // SimTrackster from CaloParticle
2823  iSim = 0;
2824 
2825  // SimTrackster with simHits connected via detid with the rechit under study
2826  //So, from all layers clusters, find the rechits that are connected with a calo particle and save/calculate the
2827  //energy of that calo particle as the sum over all rechits of the rechits energy weighted
2828  //by the caloparticle's fraction related to that rechit.
2829  const auto cpId = getCPId(simTS, iSTS, cPHandle_id, cpToSc_SimTrackstersMap, simTSs_fromCP);
2830  if (std::find(cPIndices.begin(), cPIndices.end(), cpId) == cPIndices.end())
2831  continue;
2832 
2833  CPEnergyInTS[cpId] += shared_fraction * hitEn;
2834  //Here cPOnLayer[caloparticle][layer] describe above is set.
2835  //Here for Tracksters with matched rechit the CP fraction times hit energy is added and saved .
2836  cPOnLayer[cpId].layerClusterIdToEnergyAndScore[tstId].first += shared_fraction * hitEn;
2837  sCOnLayer[cpId][iSim].layerClusterIdToEnergyAndScore[tstId].first += shared_fraction * hitEn;
2838  cPOnLayer[cpId].layerClusterIdToEnergyAndScore[tstId].second = FLT_MAX;
2839  sCOnLayer[cpId][iSim].layerClusterIdToEnergyAndScore[tstId].second = FLT_MAX;
2840  //stsInTrackster[trackster][STSids]
2841  //Connects a Trackster with all related SimTracksters.
2842  stsInTrackster[tstId].emplace_back(iSTS, FLT_MAX);
2843  //From all CaloParticles related to a layer cluster, it saves id and energy of the calo particle
2844  //that after simhit-rechit matching in layer has the maximum energy.
2845  if (shared_fraction > maxCPEnergyInTS) {
2846  //energy is used only here. cpid is saved for Tracksters
2847  maxCPEnergyInTS = CPEnergyInTS[cpId];
2848  maxCPId = cpId;
2849  }
2850  }
2851  //Keep in mind here maxCPId could be zero. So, below ask for negative not including zero to count noise.
2852  hitsToCaloParticleId[iHit] = maxCPId;
2853  }
2854 
2855  } //end of loop through rechits of the layer cluster.
2856 
2857  //Loop through all rechits to count how many of them are noise and how many are matched.
2858  //In case of matched rechit-simhit, he counts and saves the number of rechits related to the maximum energy CaloParticle.
2859  for (auto c : hitsToCaloParticleId) {
2860  if (c < 0)
2861  numberOfNoiseHitsInTS++;
2862  else
2863  occurrencesCPinTS[c]++;
2864  }
2865 
2866  //Below from all maximum energy CaloParticles, he saves the one with the largest amount
2867  //of related rechits.
2868  for (auto& c : occurrencesCPinTS) {
2869  if (c.second > maxCPNumberOfHitsInTS) {
2870  maxCPId_byNumberOfHits = c.first;
2871  maxCPNumberOfHitsInTS = c.second;
2872  }
2873  }
2874 
2875  //Find the CaloParticle that has the maximum energy shared with the Trackster under study.
2876  for (auto& c : CPEnergyInTS) {
2877  if (c.second > maxEnergySharedTSandCP) {
2878  maxCPId_byEnergy = c.first;
2879  maxEnergySharedTSandCP = c.second;
2880  }
2881  }
2882  //The energy of the CaloParticle that found to have the maximum energy shared with the Trackster under study.
2883  float totalCPEnergyFromLayerCP = 0.f;
2884  if (maxCPId_byEnergy >= 0) {
2885  totalCPEnergyFromLayerCP += cPOnLayer[maxCPId_byEnergy].energy;
2886  energyFractionOfCPinTS = maxEnergySharedTSandCP / totalCPEnergyFromLayerCP;
2887  if (tst.raw_energy() > 0.f) {
2888  energyFractionOfTSinCP = maxEnergySharedTSandCP / tst.raw_energy();
2889  }
2890  }
2891 
2892  LogDebug("HGCalValidator") << std::setw(12) << "Trackster\t" << std::setw(10) << "energy\t" << std::setw(5)
2893  << "nhits\t" << std::setw(12) << "noise hits\t" << std::setw(22)
2894  << "maxCPId_byNumberOfHits\t" << std::setw(8) << "nhitsCP\t" << std::setw(16)
2895  << "maxCPId_byEnergy\t" << std::setw(23) << "maxEnergySharedTSandCP\t" << std::setw(22)
2896  << "totalCPEnergyFromAllLayerCP\t" << std::setw(22) << "energyFractionOfTSinCP\t"
2897  << std::setw(25) << "energyFractionOfCPinTS\t" << std::endl;
2898  LogDebug("HGCalValidator") << std::setw(12) << tstId << "\t" //LogDebug("HGCalValidator")
2899  << std::setw(10) << tst.raw_energy() << "\t" << std::setw(5) << numberOfHitsInTS << "\t"
2900  << std::setw(12) << numberOfNoiseHitsInTS << "\t" << std::setw(22)
2901  << maxCPId_byNumberOfHits << "\t" << std::setw(8) << maxCPNumberOfHitsInTS << "\t"
2902  << std::setw(16) << maxCPId_byEnergy << "\t" << std::setw(23) << maxEnergySharedTSandCP
2903  << "\t" << std::setw(22) << totalCPEnergyFromLayerCP << "\t" << std::setw(22)
2904  << energyFractionOfTSinCP << "\t" << std::setw(25) << energyFractionOfCPinTS
2905  << std::endl;
2906 
2907  } // end of loop through Tracksters
2908 
2909  // Loop through Tracksters
2910  for (unsigned int tstId = 0; tstId < nTracksters; ++tstId) {
2911  const auto& tst = tracksters[tstId];
2912  if (tst.vertices().empty())
2913  continue;
2914 
2915  // find the unique SimTrackster ids contributing to the Trackster
2916  //stsInTrackster[trackster][STSids]
2917  std::sort(stsInTrackster[tstId].begin(), stsInTrackster[tstId].end());
2918  const auto last = std::unique(stsInTrackster[tstId].begin(), stsInTrackster[tstId].end());
2919  stsInTrackster[tstId].erase(last, stsInTrackster[tstId].end());
2920 
2921  if (tst.raw_energy() == 0. && !stsInTrackster[tstId].empty()) {
2922  //Loop through all SimTracksters contributing to Trackster tstId
2923  for (auto& stsPair : stsInTrackster[tstId]) {
2924  // In case of a Trackster with zero energy but related SimTracksters the score is set to 1
2925  stsPair.second = 1.;
2926  LogDebug("HGCalValidator") << "Trackster Id:\t" << tstId << "\tSimTrackster id:\t" << stsPair.first
2927  << "\tscore\t" << stsPair.second << std::endl;
2928  histograms.h_score_trackster2caloparticle[valType][count]->Fill(stsPair.second);
2929  }
2930  continue;
2931  }
2932 
2933  const auto tst_hitsAndFractions = apply_LCMultiplicity(tst, layerClusters);
2934 
2935  // Compute the correct normalization
2936  float tracksterEnergy = 0.f, invTracksterEnergyWeight = 0.f;
2937  for (const auto& haf : tst_hitsAndFractions) {
2938  float hitFr = 0.f;
2939  if (valType == 0) {
2940  hitFr = haf.second;
2941  } else {
2942  const auto lcId = getLCId(tst.vertices(), layerClusters, haf.first);
2943  const auto iLC = std::find(tst.vertices().begin(), tst.vertices().end(), lcId);
2944  hitFr = 1.f / tst.vertex_multiplicity(std::distance(std::begin(tst.vertices()), iLC));
2945  }
2946  tracksterEnergy += hitFr * hitMap.at(haf.first)->energy();
2947  invTracksterEnergyWeight += pow(hitFr * hitMap.at(haf.first)->energy(), 2);
2948  }
2949  if (invTracksterEnergyWeight)
2950  invTracksterEnergyWeight = 1.f / invTracksterEnergyWeight;
2951 
2952  for (const auto& haf : tst_hitsAndFractions) {
2953  const auto rh_detid = haf.first;
2954  unsigned int elemId = 0;
2955  float rhFraction = 0.f;
2956  if (valType == 0) {
2957  elemId = rh_detid.rawId();
2958  rhFraction = haf.second;
2959  } else {
2960  const auto lcId = getLCId(tst.vertices(), layerClusters, rh_detid);
2961  elemId = lcId;
2962  const auto iLC = std::find(tst.vertices().begin(), tst.vertices().end(), lcId);
2963  rhFraction = 1.f / tst.vertex_multiplicity(std::distance(std::begin(tst.vertices()), iLC));
2964  }
2965 
2966  bool hitWithNoSTS = false;
2967  if (detIdSimTSId_Map.find(elemId) == detIdSimTSId_Map.end())
2968  hitWithNoSTS = true;
2969  const HGCRecHit* hit = hitMap.find(rh_detid)->second;
2970  const auto hitEnergyWeight = pow(hit->energy(), 2);
2971 
2972  for (auto& stsPair : stsInTrackster[tstId]) {
2973  float cpFraction = 0.f;
2974  if (!hitWithNoSTS) {
2975  const auto& findSTSIt = std::find(
2976  detIdSimTSId_Map[elemId].begin(),
2977  detIdSimTSId_Map[elemId].end(),
2978  HGVHistoProducerAlgo::detIdInfoInCluster{
2979  stsPair.first, 0.f}); // only the first element is used for the matching (overloaded operator==)
2980  if (findSTSIt != detIdSimTSId_Map[elemId].end())
2981  cpFraction = findSTSIt->fraction;
2982  }
2983  if (stsPair.second == FLT_MAX) {
2984  stsPair.second = 0.f;
2985  }
2986  stsPair.second +=
2987  min(pow(rhFraction - cpFraction, 2), pow(rhFraction, 2)) * hitEnergyWeight * invTracksterEnergyWeight;
2988  }
2989  } // end of loop through trackster rechits
2990 
2991  //In case of a Trackster with some energy but none related CaloParticles print some info.
2992  if (stsInTrackster[tstId].empty())
2993  LogDebug("HGCalValidator") << "Trackster Id: " << tstId << "\tSimTrackster id: -1"
2994  << "\tscore: -1\n";
2995 
2996  tracksters_FakeMerge[tstId] =
2997  std::count_if(std::begin(stsInTrackster[tstId]),
2998  std::end(stsInTrackster[tstId]),
2999  [ScoreCutTStoSTSFakeMerge](const auto& obj) { return obj.second < ScoreCutTStoSTSFakeMerge; });
3000 
3001  const auto score = std::min_element(std::begin(stsInTrackster[tstId]),
3002  std::end(stsInTrackster[tstId]),
3003  [](const auto& obj1, const auto& obj2) { return obj1.second < obj2.second; });
3004  float score2 = -1;
3005  float sharedEneFrac2 = 0;
3006  for (const auto& stsPair : stsInTrackster[tstId]) {
3007  const auto iSTS = stsPair.first;
3008  const auto iScore = stsPair.second;
3009  const auto cpId = getCPId(simTSs[iSTS], iSTS, cPHandle_id, cpToSc_SimTrackstersMap, simTSs_fromCP);
3010  auto iSim = simTSs[iSTS].seedIndex();
3011  if (simTSs[iSTS].seedID() == cPHandle_id) // SimTrackster from CaloParticle
3012  iSim = 0;
3013  const auto& simOnLayer = (valType == 0) ? cPOnLayer[cpId] : sCOnLayer[cpId][iSim];
3014 
3015  float sharedeneCPallLayers = 0.;
3016  sharedeneCPallLayers += simOnLayer.layerClusterIdToEnergyAndScore.count(tstId)
3017  ? simOnLayer.layerClusterIdToEnergyAndScore.at(tstId).first
3018  : 0;
3019  if (tracksterEnergy == 0)
3020  continue;
3021  const auto sharedEneFrac = sharedeneCPallLayers / tracksterEnergy;
3022  LogDebug("HGCalValidator") << "\nTrackster id: " << tstId << " (" << tst.vertices().size() << " vertices)"
3023  << "\tSimTrackster Id: " << iSTS << " (" << simTSs[iSTS].vertices().size()
3024  << " vertices)"
3025  << " (CP id: " << cpId << ")\tscore: " << iScore
3026  << "\tsharedeneCPallLayers: " << sharedeneCPallLayers << std::endl;
3027 
3028  histograms.h_score_trackster2caloparticle[valType][count]->Fill(iScore);
3029  histograms.h_sharedenergy_trackster2caloparticle[valType][count]->Fill(sharedEneFrac);
3030  histograms.h_energy_vs_score_trackster2caloparticle[valType][count]->Fill(iScore, sharedEneFrac);
3031  if (iSTS == score->first) {
3032  histograms.h_score_trackster2bestCaloparticle[valType][count]->Fill(iScore);
3033  histograms.h_sharedenergy_trackster2bestCaloparticle[valType][count]->Fill(sharedEneFrac);
3034  histograms.h_sharedenergy_trackster2bestCaloparticle_vs_eta[valType][count]->Fill(tst.barycenter().eta(),
3035  sharedEneFrac);
3036  histograms.h_sharedenergy_trackster2bestCaloparticle_vs_phi[valType][count]->Fill(tst.barycenter().phi(),
3037  sharedEneFrac);
3038  histograms.h_energy_vs_score_trackster2bestCaloparticle[valType][count]->Fill(iScore, sharedEneFrac);
3039  } else if (score2 < 0 || iScore < score2) {
3040  score2 = iScore;
3041  sharedEneFrac2 = sharedEneFrac;
3042  }
3043  } // end of loop through SimTracksters associated to Trackster
3044  if (score2 > -1) {
3045  histograms.h_score_trackster2bestCaloparticle2[valType][count]->Fill(score2);
3046  histograms.h_sharedenergy_trackster2bestCaloparticle2[valType][count]->Fill(sharedEneFrac2);
3047  histograms.h_energy_vs_score_trackster2bestCaloparticle2[valType][count]->Fill(score2, sharedEneFrac2);
3048  }
3049  } // end of loop through Tracksters
3050 
3051  std::unordered_map<unsigned int, std::vector<float>> score3d;
3052  std::unordered_map<unsigned int, std::vector<float>> tstSharedEnergy;
3053 
3054  for (unsigned int iSTS = 0; iSTS < nSimTracksters; ++iSTS) {
3055  score3d[iSTS].resize(nTracksters);
3056  tstSharedEnergy[iSTS].resize(nTracksters);
3057  for (unsigned int j = 0; j < nTracksters; ++j) {
3058  score3d[iSTS][j] = FLT_MAX;
3059  tstSharedEnergy[iSTS][j] = 0.f;
3060  }
3061  }
3062 
3063  // Fill the plots to compute the different metrics linked to
3064  // gen-level, namely efficiency, purity and duplicate. In this loop should restrict
3065  // only to the selected caloParaticles.
3066  for (unsigned int iSTS = 0; iSTS < nSimTracksters; ++iSTS) {
3067  const auto& sts = simTSs[iSTS];
3068  const auto& cpId = getCPId(sts, iSTS, cPHandle_id, cpToSc_SimTrackstersMap, simTSs_fromCP);
3069  if (valType == 0 && std::find(cPSelectedIndices.begin(), cPSelectedIndices.end(), cpId) == cPSelectedIndices.end())
3070  continue;
3071 
3072  const auto& hafLC = apply_LCMultiplicity(sts, layerClusters);
3073  float SimEnergy_LC = 0.f;
3074  for (const auto& haf : hafLC) {
3075  const auto lcId = getLCId(sts.vertices(), layerClusters, haf.first);
3076  const auto iLC = std::find(sts.vertices().begin(), sts.vertices().end(), lcId);
3077  SimEnergy_LC +=
3078  hitMap.at(haf.first)->energy() / sts.vertex_multiplicity(std::distance(std::begin(sts.vertices()), iLC));
3079  }
3080 
3081  auto iSim = sts.seedIndex();
3082  if (sts.seedID() == cPHandle_id) // SimTrackster from CaloParticle
3083  iSim = 0;
3084  auto& simOnLayer = (valType == 0) ? cPOnLayer[cpId] : sCOnLayer[cpId][iSim];
3085 
3086  // Keep the Trackster ids that are related to
3087  // SimTrackster under study for the final filling of the score
3088  std::set<unsigned int> stsId_tstId_related;
3089  auto& score3d_iSTS = score3d[iSTS];
3090 
3091  float SimEnergy = 0.f;
3092  float SimEnergyWeight = 0.f, hitsEnergyWeight = 0.f;
3093  //for (unsigned int layerId = 0; layerId < 1/*layers * 2*/; ++layerId) {
3094  const auto SimNumberOfHits = simOnLayer.hits_and_fractions.size();
3095  if (SimNumberOfHits == 0)
3096  continue;
3097  SimEnergy += simOnLayer.energy;
3098  int tstWithMaxEnergyInCP = -1;
3099  //This is the maximum energy related to Trackster per layer.
3100  float maxEnergyTSperlayerinSim = 0.f;
3101  float SimEnergyFractionInTSperlayer = 0.f;
3102  //Remember and not confused by name. layerClusterIdToEnergyAndScore contains the Trackster id.
3103  for (const auto& tst : simOnLayer.layerClusterIdToEnergyAndScore) {
3104  if (tst.second.first > maxEnergyTSperlayerinSim) {
3105  maxEnergyTSperlayerinSim = tst.second.first;
3106  tstWithMaxEnergyInCP = tst.first;
3107  }
3108  }
3109  if (SimEnergy > 0.f)
3110  SimEnergyFractionInTSperlayer = maxEnergyTSperlayerinSim / SimEnergy;
3111 
3112  LogDebug("HGCalValidator") << std::setw(12) << "caloparticle\t" << std::setw(15) << "cp total energy\t"
3113  << std::setw(15) << "cpEnergyOnLayer\t" << std::setw(14) << "CPNhitsOnLayer\t"
3114  << std::setw(18) << "tstWithMaxEnergyInCP\t" << std::setw(15) << "maxEnergyTSinCP\t"
3115  << std::setw(20) << "CPEnergyFractionInTS"
3116  << "\n";
3117  LogDebug("HGCalValidator") << std::setw(12) << cpId << "\t" << std::setw(15) << sts.raw_energy() << "\t"
3118  << std::setw(15) << SimEnergy << "\t" << std::setw(14) << SimNumberOfHits << "\t"
3119  << std::setw(18) << tstWithMaxEnergyInCP << "\t" << std::setw(15)
3120  << maxEnergyTSperlayerinSim << "\t" << std::setw(20) << SimEnergyFractionInTSperlayer
3121  << "\n";
3122 
3123  for (const auto& haf : ((valType == 0) ? simOnLayer.hits_and_fractions : hafLC)) {
3124  const auto& hitDetId = haf.first;
3125  // Compute the correct normalization
3126  // Need to loop on the simOnLayer data structure since this is the
3127  // only one that has the compressed information for multiple usage
3128  // of the same DetId by different SimClusters by a single CaloParticle.
3129  SimEnergyWeight += pow(haf.second * hitMap.at(hitDetId)->energy(), 2);
3130 
3131  const auto lcId = getLCId(sts.vertices(), layerClusters, hitDetId);
3132  float cpFraction = 0.f;
3133  if (valType == 0) {
3134  cpFraction = haf.second;
3135  } else {
3136  const auto iLC = std::find(sts.vertices().begin(), sts.vertices().end(), lcId);
3137  cpFraction = 1.f / sts.vertex_multiplicity(std::distance(std::begin(sts.vertices()), iLC));
3138  }
3139  if (cpFraction == 0.f)
3140  continue; // hopefully this should never happen
3141 
3142  bool hitWithNoTS = false;
3143  if (detIdToTracksterId_Map.find(hitDetId) == detIdToTracksterId_Map.end())
3144  hitWithNoTS = true;
3145  const HGCRecHit* hit = hitMap.find(hitDetId)->second;
3146  const auto hitEnergyWeight = pow(hit->energy(), 2);
3147  hitsEnergyWeight += pow(cpFraction, 2) * hitEnergyWeight;
3148 
3149  for (auto& tsPair : simOnLayer.layerClusterIdToEnergyAndScore) {
3150  const auto tstId = tsPair.first;
3151  stsId_tstId_related.insert(tstId);
3152 
3153  float tstFraction = 0.f;
3154  if (!hitWithNoTS) {
3155  const auto findTSIt =
3156  std::find(detIdToTracksterId_Map[hitDetId].begin(),
3157  detIdToTracksterId_Map[hitDetId].end(),
3158  HGVHistoProducerAlgo::detIdInfoInTrackster{
3159  tstId, 0, 0.f}); // only the first element is used for the matching (overloaded operator==)
3160  if (findTSIt != detIdToTracksterId_Map[hitDetId].end()) {
3161  if (valType == 0) {
3162  tstFraction = findTSIt->fraction;
3163  } else {
3164  const auto iLC = std::find(
3165  tracksters[tstId].vertices().begin(), tracksters[tstId].vertices().end(), findTSIt->clusterId);
3166  if (iLC != tracksters[tstId].vertices().end()) {
3167  tstFraction = 1.f / tracksters[tstId].vertex_multiplicity(
3168  std::distance(std::begin(tracksters[tstId].vertices()), iLC));
3169  }
3170  }
3171  }
3172  }
3173  // Here do not divide as before by the trackster energy weight. Should sum first
3174  // over all layers and divide with the total CP energy over all layers.
3175  if (tsPair.second.second == FLT_MAX) {
3176  tsPair.second.second = 0.f;
3177  }
3178  tsPair.second.second += min(pow(tstFraction - cpFraction, 2), pow(cpFraction, 2)) * hitEnergyWeight;
3179 
3180  LogDebug("HGCalValidator") << "\nTracksterId:\t" << tstId << "\tSimTracksterId:\t" << iSTS << "\tcpId:\t"
3181  << cpId << "\ttstfraction, cpfraction:\t" << tstFraction << ", " << cpFraction
3182  << "\thitEnergyWeight:\t" << hitEnergyWeight << "\tadded delta:\t"
3183  << pow((tstFraction - cpFraction), 2) * hitEnergyWeight
3184  << "\tcurrent Sim-score numerator:\t" << tsPair.second.second
3185  << "\tshared Sim energy:\t" << tsPair.second.first << '\n';
3186  }
3187  } // end of loop through SimCluster SimHits on current layer
3188 
3189  if (simOnLayer.layerClusterIdToEnergyAndScore.empty())
3190  LogDebug("HGCalValidator") << "CP Id:\t" << cpId << "\tTS id:\t-1"
3191  << " Sub score in \t -1\n";
3192 
3193  for (const auto& tsPair : simOnLayer.layerClusterIdToEnergyAndScore) {
3194  const auto tstId = tsPair.first;
3195  // 3D score here without the denominator at this point
3196  if (score3d_iSTS[tstId] == FLT_MAX) {
3197  score3d_iSTS[tstId] = 0.f;
3198  }
3199  score3d_iSTS[tstId] += tsPair.second.second;
3200  tstSharedEnergy[iSTS][tstId] += tsPair.second.first;
3201  }
3202  //} // end of loop through layers
3203 
3204  const auto scoreDenom = (valType == 0) ? SimEnergyWeight : hitsEnergyWeight;
3205  const auto energyDenom = (valType == 0) ? SimEnergy : SimEnergy_LC;
3206 
3207  const auto sts_eta = sts.barycenter().eta();
3208  const auto sts_phi = sts.barycenter().phi();
3209  const auto sts_en = sts.raw_energy();
3210  const auto sts_pt = sts.raw_pt();
3211  histograms.h_denom_caloparticle_eta[valType][count]->Fill(sts_eta);
3212  histograms.h_denom_caloparticle_phi[valType][count]->Fill(sts_phi);
3213  histograms.h_denom_caloparticle_en[valType][count]->Fill(sts_en);
3214  histograms.h_denom_caloparticle_pt[valType][count]->Fill(sts_pt);
3215 
3216  //Loop through related Tracksters here
3217  // In case the threshold to associate a CaloParticle to a Trackster is
3218  // below 50%, there could be cases in which the CP is linked to more than
3219  // one tracksters, leading to efficiencies >1. This boolean is used to
3220  // avoid "over counting".
3221  bool sts_considered_efficient = false;
3222  bool sts_considered_pure = false;
3223  for (const auto tstId : stsId_tstId_related) {
3224  // Now time for the denominator
3225  score3d_iSTS[tstId] /= scoreDenom;
3226  const auto tstSharedEnergyFrac = tstSharedEnergy[iSTS][tstId] / energyDenom;
3227  LogDebug("HGCalValidator") << "STS id: " << iSTS << "\t(CP id: " << cpId << ")\tTS id: " << tstId
3228  << "\nSimEnergy: " << energyDenom << "\tSimEnergyWeight: " << SimEnergyWeight
3229  << "\tTrackste energy: " << tracksters[tstId].raw_energy()
3230  << "\nscore: " << score3d_iSTS[tstId]
3231  << "\tshared energy: " << tstSharedEnergy[iSTS][tstId]
3232  << "\tshared energy fraction: " << tstSharedEnergyFrac << "\n";
3233 
3234  histograms.h_score_caloparticle2trackster[valType][count]->Fill(score3d_iSTS[tstId]);
3235  histograms.h_sharedenergy_caloparticle2trackster[valType][count]->Fill(tstSharedEnergyFrac);
3236  histograms.h_energy_vs_score_caloparticle2trackster[valType][count]->Fill(score3d_iSTS[tstId],
3237  tstSharedEnergyFrac);
3238  // Fill the numerator for the efficiency calculation. The efficiency is computed by considering the energy shared between a Trackster and a _corresponding_ caloParticle. The threshold is configurable via python.
3239  if (!sts_considered_efficient && (tstSharedEnergyFrac >= minTSTSharedEneFracEfficiency_)) {
3240  sts_considered_efficient = true;
3241  histograms.h_numEff_caloparticle_eta[valType][count]->Fill(sts_eta);
3242  histograms.h_numEff_caloparticle_phi[valType][count]->Fill(sts_phi);
3243  histograms.h_numEff_caloparticle_en[valType][count]->Fill(sts_en);
3244  histograms.h_numEff_caloparticle_pt[valType][count]->Fill(sts_pt);
3245  }
3246 
3247  if (score3d_iSTS[tstId] < ScoreCutSTStoTSPurDup) {
3248  if (tracksters_PurityDuplicate[tstId] < 1)
3249  tracksters_PurityDuplicate[tstId]++; // for Purity
3250  if (sts_considered_pure)
3251  tracksters_PurityDuplicate[tstId]++; // for Duplicate
3252  sts_considered_pure = true;
3253  }
3254  } // end of loop through Tracksters related to SimTrackster
3255 
3256  const auto best = std::min_element(std::begin(score3d_iSTS), std::end(score3d_iSTS));
3257  if (best != score3d_iSTS.end()) {
3258  const auto bestTstId = std::distance(std::begin(score3d_iSTS), best);
3259  const auto bestTstSharedEnergyFrac = tstSharedEnergy[iSTS][bestTstId] / energyDenom;
3260  histograms.h_scorePur_caloparticle2trackster[valType][count]->Fill(*best);
3261  histograms.h_sharedenergy_caloparticle2trackster_assoc[valType][count]->Fill(bestTstSharedEnergyFrac);
3262  histograms.h_sharedenergy_caloparticle2trackster_assoc_vs_eta[valType][count]->Fill(sts_eta,
3263  bestTstSharedEnergyFrac);
3264  histograms.h_sharedenergy_caloparticle2trackster_assoc_vs_phi[valType][count]->Fill(sts_phi,
3265  bestTstSharedEnergyFrac);
3266  histograms.h_energy_vs_score_caloparticle2bestTrackster[valType][count]->Fill(*best, bestTstSharedEnergyFrac);
3267  LogDebug("HGCalValidator") << count << " " << sts_eta << " " << sts_phi << " "
3268  << tracksters[bestTstId].raw_energy() << " " << sts.raw_energy() << " "
3269  << bestTstSharedEnergyFrac << "\n";
3270 
3271  if (score3d_iSTS.size() > 1) {
3272  auto best2 = (best == score3d_iSTS.begin()) ? std::next(best, 1) : score3d_iSTS.begin();
3273  for (auto tstId = score3d_iSTS.begin(); tstId != score3d_iSTS.end() && tstId != best; tstId++)
3274  if (*tstId < *best2)
3275  best2 = tstId;
3276  const auto best2TstId = std::distance(std::begin(score3d_iSTS), best2);
3277  const auto best2TstSharedEnergyFrac = tstSharedEnergy[iSTS][best2TstId] / energyDenom;
3278  histograms.h_scoreDupl_caloparticle2trackster[valType][count]->Fill(*best2);
3279  histograms.h_sharedenergy_caloparticle2trackster_assoc2[valType][count]->Fill(best2TstSharedEnergyFrac);
3280  histograms.h_energy_vs_score_caloparticle2bestTrackster2[valType][count]->Fill(*best2,
3281  best2TstSharedEnergyFrac);
3282  }
3283  }
3284  } // end of loop through SimTracksters
3285 
3286  // Fill the plots to compute the different metrics linked to
3287  // reco-level, namely fake-rate an merge-rate. Should *not*
3288  // restrict only to the selected caloParaticles.
3289  for (unsigned int tstId = 0; tstId < nTracksters; ++tstId) {
3290  const auto& tst = tracksters[tstId];
3291  if (tst.vertices().empty())
3292  continue;
3293  const auto iTS_eta = tst.barycenter().eta();
3294  const auto iTS_phi = tst.barycenter().phi();
3295  const auto iTS_en = tst.raw_energy();
3296  const auto iTS_pt = tst.raw_pt();
3297  histograms.h_denom_trackster_eta[valType][count]->Fill(iTS_eta);
3298  histograms.h_denom_trackster_phi[valType][count]->Fill(iTS_phi);
3299  histograms.h_denom_trackster_en[valType][count]->Fill(iTS_en);
3300  histograms.h_denom_trackster_pt[valType][count]->Fill(iTS_pt);
3301 
3302  if (tracksters_PurityDuplicate[tstId] > 0) {
3303  histograms.h_num_caloparticle_eta[valType][count]->Fill(iTS_eta);
3304  histograms.h_num_caloparticle_phi[valType][count]->Fill(iTS_phi);
3305  histograms.h_num_caloparticle_en[valType][count]->Fill(iTS_en);
3306  histograms.h_num_caloparticle_pt[valType][count]->Fill(iTS_pt);
3307 
3308  if (tracksters_PurityDuplicate[tstId] > 1) {
3309  histograms.h_numDup_trackster_eta[valType][count]->Fill(iTS_eta);
3310  histograms.h_numDup_trackster_phi[valType][count]->Fill(iTS_phi);
3311  histograms.h_numDup_trackster_en[valType][count]->Fill(iTS_en);
3312  histograms.h_numDup_trackster_pt[valType][count]->Fill(iTS_pt);
3313  }
3314  }
3315 
3316  if (tracksters_FakeMerge[tstId] > 0) {
3317  histograms.h_num_trackster_eta[valType][count]->Fill(iTS_eta);
3318  histograms.h_num_trackster_phi[valType][count]->Fill(iTS_phi);
3319  histograms.h_num_trackster_en[valType][count]->Fill(iTS_en);
3320  histograms.h_num_trackster_pt[valType][count]->Fill(iTS_pt);
3321 
3322  if (tracksters_FakeMerge[tstId] > 1) {
3323  histograms.h_numMerge_trackster_eta[valType][count]->Fill(iTS_eta);
3324  histograms.h_numMerge_trackster_phi[valType][count]->Fill(iTS_phi);
3325  histograms.h_numMerge_trackster_en[valType][count]->Fill(iTS_en);
3326  histograms.h_numMerge_trackster_pt[valType][count]->Fill(iTS_pt);
3327  }
3328  }
3329  } // End loop over Tracksters
3330 }
3331 
3332 void HGVHistoProducerAlgo::fill_trackster_histos(
3333  const Histograms& histograms,
3334  const int count,
3335  const ticl::TracksterCollection& tracksters,
3336  const reco::CaloClusterCollection& layerClusters,
3337  const ticl::TracksterCollection& simTSs,
3338  const ticl::TracksterCollection& simTSs_fromCP,
3339  const std::map<unsigned int, std::vector<unsigned int>>& cpToSc_SimTrackstersMap,
3340  std::vector<SimCluster> const& sC,
3341  const edm::ProductID& cPHandle_id,
3342  std::vector<CaloParticle> const& cP,
3343  std::vector<size_t> const& cPIndices,
3344  std::vector<size_t> const& cPSelectedIndices,
3345  std::unordered_map<DetId, const HGCRecHit*> const& hitMap,
3346  unsigned int layers) const {
3347  //Each event to be treated as two events:
3348  //an event in +ve endcap, plus another event in -ve endcap.
3349 
3350  //To keep track of total num of Tracksters
3351  int totNTstZm = 0; //-z
3352  int totNTstZp = 0; //+z
3353  //To count the number of Tracksters with 3 contiguous layers per event.
3354  int totNContTstZp = 0; //+z
3355  int totNContTstZm = 0; //-z
3356  //For the number of Tracksters without 3 contiguous layers per event.
3357  int totNNotContTstZp = 0; //+z
3358  int totNNotContTstZm = 0; //-z
3359  // Check below the score of cont and non cont Tracksters
3360  std::vector<bool> contTracksters;
3361  contTracksters.clear();
3362 
3363  //[tstId]-> vector of 2d layer clusters size
3364  std::unordered_map<unsigned int, std::vector<unsigned int>> multiplicity;
3365  //[tstId]-> [layer][cluster size]
3366  std::unordered_map<unsigned int, std::vector<unsigned int>> multiplicity_vs_layer;
3367  //We will need for the scale text option
3368  // unsigned int totalLcInTsts = 0;
3369  // for (unsigned int tstId = 0; tstId < nTracksters; ++tstId) {
3370  // totalLcInTsts = totalLcInTsts + tracksters[tstId].vertices().size();
3371  // }
3372 
3373  const auto nTracksters = tracksters.size();
3374  // loop through Tracksters
3375  for (unsigned int tstId = 0; tstId < nTracksters; ++tstId) {
3376  const auto& tst = tracksters[tstId];
3377  if (tst.vertices().empty())
3378  continue;
3379 
3380  if (tst.barycenter().z() < 0.)
3381  totNTstZm++;
3382  else if (tst.barycenter().z() > 0.)
3383  totNTstZp++;
3384 
3385  //Total number of layer clusters in Trackster
3386  int tnLcInTst = 0;
3387 
3388  //To keep track of total num of layer clusters per Trackster
3389  //tnLcInTstperlaypz[layerid], tnLcInTstperlaymz[layerid]
3390  std::vector<int> tnLcInTstperlay(1000, 0); //+z
3391 
3392  //For the layers the Trackster expands to. Will use a set because there would be many
3393  //duplicates and then go back to vector for random access, since they say it is faster.
3394  std::set<unsigned int> trackster_layers;
3395 
3396  bool tracksterInZplus = false;
3397  bool tracksterInZminus = false;
3398 
3399  //Loop through layer clusters
3400  for (const auto lcId : tst.vertices()) {
3401  //take the hits and their fraction of the specific layer cluster.
3402  const auto& hits_and_fractions = layerClusters[lcId].hitsAndFractions();
3403 
3404  //For the multiplicity of the 2d layer clusters in Tracksters
3405  multiplicity[tstId].emplace_back(hits_and_fractions.size());
3406 
3407  const auto firstHitDetId = hits_and_fractions[0].first;
3408  //The layer that the layer cluster belongs to
3409  const auto layerid = recHitTools_->getLayerWithOffset(firstHitDetId) +
3410  layers * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
3411  trackster_layers.insert(layerid);
3412  multiplicity_vs_layer[tstId].emplace_back(layerid);
3413 
3414  tnLcInTstperlay[layerid]++;
3415  tnLcInTst++;
3416 
3417  if (recHitTools_->zside(firstHitDetId) > 0.)
3418  tracksterInZplus = true;
3419  else if (recHitTools_->zside(firstHitDetId) < 0.)
3420  tracksterInZminus = true;
3421  } // end of loop through layerClusters
3422 
3423  // Per layer : Loop 0->99
3424  for (unsigned ilayer = 0; ilayer < layers * 2; ++ilayer) {
3425  if (histograms.h_clusternum_in_trackster_perlayer[count].count(ilayer) && tnLcInTstperlay[ilayer] != 0) {
3426  histograms.h_clusternum_in_trackster_perlayer[count].at(ilayer)->Fill((float)tnLcInTstperlay[ilayer]);
3427  }
3428  // For the profile now of 2d layer cluster in Tracksters vs layer number.
3429  if (tnLcInTstperlay[ilayer] != 0) {
3430  histograms.h_clusternum_in_trackster_vs_layer[count]->Fill((float)ilayer, (float)tnLcInTstperlay[ilayer]);
3431  }
3432  } // end of loop over layers
3433 
3434  // Looking for Tracksters with 3 contiguous layers per event.
3435  std::vector<int> trackster_layers_vec(trackster_layers.begin(), trackster_layers.end());
3436  // Check also for non contiguous Tracksters
3437  bool contiTrackster = false;
3438  // Start from 1 and go up to size - 1 element.
3439  if (trackster_layers_vec.size() >= 3) {
3440  for (unsigned int iLayer = 1; iLayer < trackster_layers_vec.size() - 1; ++iLayer) {
3441  if ((trackster_layers_vec[iLayer - 1] + 1 == trackster_layers_vec[iLayer]) &&
3442  (trackster_layers_vec[iLayer + 1] - 1 == trackster_layers_vec[iLayer])) {
3443  // Trackster with 3 contiguous layers per event
3444  if (tracksterInZplus)
3445  totNContTstZp++;
3446  else if (tracksterInZminus)
3447  totNContTstZm++;
3448 
3449  contiTrackster = true;
3450  break;
3451  }
3452  }
3453  }
3454  // Count non contiguous Tracksters
3455  if (!contiTrackster) {
3456  if (tracksterInZplus)
3457  totNNotContTstZp++;
3458  else if (tracksterInZminus)
3459  totNNotContTstZm++;
3460  }
3461 
3462  // Save for the score
3463  contTracksters.push_back(contiTrackster);
3464 
3465  histograms.h_clusternum_in_trackster[count]->Fill(tnLcInTst);
3466 
3467  for (unsigned int lc = 0; lc < multiplicity[tstId].size(); ++lc) {
3468  //multiplicity of the current LC
3469  float mlp = std::count(std::begin(multiplicity[tstId]), std::end(multiplicity[tstId]), multiplicity[tstId][lc]);
3470  //LogDebug("HGCalValidator") << "mlp %" << (100. * mlp)/ ((float) nLayerClusters) << std::endl;
3471  // histograms.h_multiplicityOfLCinTST[count]->Fill( mlp , multiplicity[tstId][lc] , 100. / (float) totalLcInTsts );
3472  histograms.h_multiplicityOfLCinTST[count]->Fill(mlp, multiplicity[tstId][lc]);
3473  //When plotting with the text option we want the entries to be the same
3474  //as the % of the current cell over the whole number of layerClusters. For this we need an extra histo.
3475  histograms.h_multiplicity_numberOfEventsHistogram[count]->Fill(mlp);
3476  //For the cluster multiplicity vs layer
3477  //First with the -z endcap (V10:0->49)
3478  if (multiplicity_vs_layer[tstId][lc] < layers) {
3479  histograms.h_multiplicityOfLCinTST_vs_layercluster_zminus[count]->Fill(mlp, multiplicity_vs_layer[tstId][lc]);
3480  histograms.h_multiplicity_zminus_numberOfEventsHistogram[count]->Fill(mlp);
3481  } else { //Then for the +z (V10:50->99)
3482  histograms.h_multiplicityOfLCinTST_vs_layercluster_zplus[count]->Fill(
3483  mlp, multiplicity_vs_layer[tstId][lc] - layers);
3484  histograms.h_multiplicity_zplus_numberOfEventsHistogram[count]->Fill(mlp);
3485  }
3486  //For the cluster multiplicity vs cluster energy
3487  histograms.h_multiplicityOfLCinTST_vs_layerclusterenergy[count]->Fill(mlp,
3488  layerClusters[tst.vertices(lc)].energy());
3489  }
3490 
3491  if (!trackster_layers.empty()) {
3492  histograms.h_trackster_x[count]->Fill(tst.barycenter().x());
3493  histograms.h_trackster_y[count]->Fill(tst.barycenter().y());
3494  histograms.h_trackster_z[count]->Fill(tst.barycenter().z());
3495  histograms.h_trackster_eta[count]->Fill(tst.barycenter().eta());
3496  histograms.h_trackster_phi[count]->Fill(tst.barycenter().phi());
3497 
3498  histograms.h_trackster_firstlayer[count]->Fill((float)*trackster_layers.begin());
3499  histograms.h_trackster_lastlayer[count]->Fill((float)*trackster_layers.rbegin());
3500  histograms.h_trackster_layersnum[count]->Fill((float)trackster_layers.size());
3501 
3502  histograms.h_trackster_pt[count]->Fill(tst.raw_pt());
3503  histograms.h_trackster_energy[count]->Fill(tst.raw_energy());
3504  }
3505 
3506  } //end of loop through Tracksters
3507 
3508  histograms.h_tracksternum[count]->Fill(totNTstZm + totNTstZp);
3509  histograms.h_conttracksternum[count]->Fill(totNContTstZp + totNContTstZm);
3510  histograms.h_nonconttracksternum[count]->Fill(totNNotContTstZp + totNNotContTstZm);
3511 
3512  // CaloParticle
3513  tracksters_to_SimTracksters(histograms,
3514  count,
3515  tracksters,
3516  layerClusters,
3517  simTSs_fromCP,
3518  Linking,
3519  simTSs_fromCP,
3520  cpToSc_SimTrackstersMap,
3521  sC,
3522  cPHandle_id,
3523  cP,
3524  cPIndices,
3525  cPSelectedIndices,
3526  hitMap,
3527  layers);
3528 
3529  // Pattern recognition
3530  tracksters_to_SimTracksters(histograms,
3531  count,
3532  tracksters,
3533  layerClusters,
3534  simTSs,
3535  PatternRecognition,
3536  simTSs_fromCP,
3537  cpToSc_SimTrackstersMap,
3538  sC,
3539  cPHandle_id,
3540  cP,
3541  cPIndices,
3542  cPSelectedIndices,
3543  hitMap,
3544  layers);
3545 }
3546 
3547 double HGVHistoProducerAlgo::distance2(const double x1,
3548  const double y1,
3549  const double x2,
3550  const double y2) const { //distance squared
3551  const double dx = x1 - x2;
3552  const double dy = y1 - y2;
3553  return (dx * dx + dy * dy);
3554 } //distance squaredq
3555 double HGVHistoProducerAlgo::distance(const double x1,
3556  const double y1,
3557  const double x2,
3558  const double y2) const { //2-d distance on the layer (x-y)
3559  return std::sqrt(distance2(x1, y1, x2, y2));
3560 }
3561 
3562 void HGVHistoProducerAlgo::setRecHitTools(std::shared_ptr<hgcal::RecHitTools> recHitTools) {
3563  recHitTools_ = recHitTools;
3564 }
3565 
3566 DetId HGVHistoProducerAlgo::findmaxhit(const reco::CaloCluster& cluster,
3567  std::unordered_map<DetId, const HGCRecHit*> const& hitMap) const {
3568  const auto& hits_and_fractions = cluster.hitsAndFractions();
3569 
3570  DetId themaxid;
3571  double maxene = 0.;
3572  for (std::vector<std::pair<DetId, float>>::const_iterator it_haf = hits_and_fractions.begin();
3573  it_haf != hits_and_fractions.end();
3574  ++it_haf) {
3575  const DetId rh_detid = it_haf->first;
3576  const auto hitEn = hitMap.find(rh_detid)->second->energy();
3577 
3578  if (maxene < hitEn) {
3579  maxene = hitEn;
3580  themaxid = rh_detid;
3581  }
3582  }
3583 
3584  return themaxid;
3585 }
3586 
3587 double HGVHistoProducerAlgo::getEta(double eta) const {
3588  if (useFabsEta_)
3589  return fabs(eta);
3590  else
3591  return eta;
3592 }
Calorimetry_cff.thickness
thickness
Definition: Calorimetry_cff.py:115
Phase2L1GMT::maxPhi_
const float maxPhi_
Definition: Constants.h:69
CaloParticle::g4Tracks
const std::vector< SimTrack > & g4Tracks() const
Definition: CaloParticle.h:74
reco::CaloCluster::hitsAndFractions
const std::vector< std::pair< DetId, float > > & hitsAndFractions() const
Definition: CaloCluster.h:210
mps_fire.i
i
Definition: mps_fire.py:429
ScoreCutSCtoLC_
const double ScoreCutSCtoLC_
Definition: HGVHistoProducerAlgo.cc:19
CaloParticle::simClusters
const SimClusterRefVector & simClusters() const
Definition: CaloParticle.h:72
heavyIonCSV_trainingSettings.idx
idx
Definition: heavyIonCSV_trainingSettings.py:5
HGVHistoProducerAlgo::bookTracksterHistos
void bookTracksterHistos(DQMStore::IBooker &ibook, Histograms &histograms, unsigned int layers)
Definition: HGVHistoProducerAlgo.cc:972
offlineSlimmedPrimaryVertices_cfi.score
score
Definition: offlineSlimmedPrimaryVertices_cfi.py:6
HGVHistoProducerAlgo::setRecHitTools
void setRecHitTools(std::shared_ptr< hgcal::RecHitTools > recHitTools)
Definition: HGVHistoProducerAlgo.cc:3562
testProducerWithPsetDescEmpty_cfi.y2
y2
Definition: testProducerWithPsetDescEmpty_cfi.py:30
HLT_2022v15_cff.distance
distance
Definition: HLT_2022v15_cff.py:6339
muonDTDigis_cfi.pset
pset
Definition: muonDTDigis_cfi.py:27
ScoreCutLCtoSC_
const double ScoreCutLCtoSC_
Definition: HGVHistoProducerAlgo.cc:18
cond::impl::to_string
std::string to_string(const V &value)
Definition: OMSAccess.h:71
hltEgammaHGCALIDVarsL1Seeded_cfi.layerClusters
layerClusters
Definition: hltEgammaHGCALIDVarsL1Seeded_cfi.py:5
weight
Definition: weight.py:1
HGVHistoProducerAlgo::fill_info_histos
void fill_info_histos(const Histograms &histograms, unsigned int layers) const
Definition: HGVHistoProducerAlgo.cc:1427
HGVHistoProducerAlgo::bookClusterHistos_LCtoCP_association
void bookClusterHistos_LCtoCP_association(DQMStore::IBooker &ibook, Histograms &histograms, unsigned int layers)
Definition: HGVHistoProducerAlgo.cc:697
EgammaValidation_cff.pdgid
pdgid
Definition: EgammaValidation_cff.py:29
ecaldqm::zside
int zside(DetId const &)
Definition: EcalDQMCommonUtils.cc:189
HLT_2022v15_cff.fraction
fraction
Definition: HLT_2022v15_cff.py:35229
spr::find
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
GetRecoTauVFromDQM_MC_cff.next
next
Definition: GetRecoTauVFromDQM_MC_cff.py:31
DetId::det
constexpr Detector det() const
get the detector field from this detid
Definition: DetId.h:46
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
gpuClustering::pixelStatus::mask
constexpr uint32_t mask
Definition: gpuClustering.h:24
edm::AssociationMap::find
const_iterator find(const key_type &k) const
find element with specified reference key
Definition: AssociationMap.h:175
HGVHistoProducerAlgo::bookInfo
void bookInfo(DQMStore::IBooker &ibook, Histograms &histograms)
Definition: HGVHistoProducerAlgo.cc:205
edm::AssociationMap::end
const_iterator end() const
last iterator over the map (read only)
Definition: AssociationMap.h:173
AlignmentTracksFromVertexSelector_cfi.vertices
vertices
Definition: AlignmentTracksFromVertexSelector_cfi.py:5
CaloParticle::energy
float energy() const
Energy. Note this is taken from the first SimTrack only.
Definition: CaloParticle.h:98
HGVHistoProducerAlgo::bookSimClusterAssociationHistos
void bookSimClusterAssociationHistos(DQMStore::IBooker &ibook, Histograms &histograms, unsigned int layers, std::vector< int > thicknesses)
Definition: HGVHistoProducerAlgo.cc:343
HGVHistoProducerAlgo::findmaxhit
DetId findmaxhit(const reco::CaloCluster &cluster, std::unordered_map< DetId, const HGCRecHit *> const &) const
Definition: HGVHistoProducerAlgo.cc:3566
HGVHistoProducerAlgo::fill_caloparticle_histos
void fill_caloparticle_histos(const Histograms &histograms, int pdgid, const CaloParticle &caloparticle, std::vector< SimVertex > const &simVertices, unsigned int layers, std::unordered_map< DetId, const HGCRecHit *> const &) const
Definition: HGVHistoProducerAlgo.cc:1440
HGVHistoProducerAlgo::bookSimClusterHistos
void bookSimClusterHistos(DQMStore::IBooker &ibook, Histograms &histograms, unsigned int layers, std::vector< int > thicknesses)
Definition: HGVHistoProducerAlgo.cc:291
HGVHistoProducerAlgo::fill_cluster_histos
void fill_cluster_histos(const Histograms &histograms, const int count, const reco::CaloCluster &cluster) const
Definition: HGVHistoProducerAlgo.cc:1732
dqm::implementation::IBooker::bookProfile
MonitorElement * bookProfile(TString const &name, TString const &title, int nchX, double lowX, double highX, int, double lowY, double highY, char const *option="s", FUNC onbooking=NOOP())
Definition: DQMStore.h:399
testProducerWithPsetDescEmpty_cfi.x2
x2
Definition: testProducerWithPsetDescEmpty_cfi.py:28
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
CaloParticle::pt
float pt() const
Transverse momentum. Note this is taken from the first SimTrack only.
Definition: CaloParticle.h:130
tier0.unique
def unique(seq, keepstr=True)
Definition: tier0.py:24
HGVHistoProducerAlgo::distance2
double distance2(const double x1, const double y1, const double x2, const double y2) const
Definition: HGVHistoProducerAlgo.cc:3547
jetsAK4_CHS_cff.sort
sort
Definition: jetsAK4_CHS_cff.py:39
HGVHistoProducerAlgo::layerClusters_to_SimClusters
void layerClusters_to_SimClusters(const Histograms &histograms, const int count, edm::Handle< reco::CaloClusterCollection > clusterHandle, const reco::CaloClusterCollection &clusters, edm::Handle< std::vector< SimCluster >> simClusterHandle, std::vector< SimCluster > const &simClusters, std::vector< size_t > const &sCIndices, const std::vector< float > &mask, std::unordered_map< DetId, const HGCRecHit *> const &, unsigned int layers, const hgcal::RecoToSimCollectionWithSimClusters &recSimColl, const hgcal::SimToRecoCollectionWithSimClusters &simRecColl) const
Definition: HGVHistoProducerAlgo.cc:1997
reco::CaloClusterCollection
std::vector< CaloCluster > CaloClusterCollection
collection of CaloCluster objects
Definition: CaloClusterFwd.h:19
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
HGVHistoProducerAlgo::bookCaloParticleHistos
void bookCaloParticleHistos(DQMStore::IBooker &ibook, Histograms &histograms, int pdgid, unsigned int layers)
Definition: HGVHistoProducerAlgo.cc:214
HGVHistoProducerAlgo::bookTracksterSTSHistos
void bookTracksterSTSHistos(DQMStore::IBooker &ibook, Histograms &histograms, const validationType valType)
Definition: HGVHistoProducerAlgo.cc:1114
ScoreCutTStoSTSFakeMerge_
const double ScoreCutTStoSTSFakeMerge_[]
Definition: HGVHistoProducerAlgo.cc:20
ScoreCutLCtoCP_
const double ScoreCutLCtoCP_
Definition: HGVHistoProducerAlgo.cc:16
CaloParticle::phi
float phi() const
Momentum azimuthal angle. Note this is taken from the first SimTrack only.
Definition: CaloParticle.h:134
HGVHistoProducerAlgo::fill_generic_cluster_histos
void fill_generic_cluster_histos(const Histograms &histograms, const int count, edm::Handle< reco::CaloClusterCollection > clusterHandle, const reco::CaloClusterCollection &clusters, const Density &densities, edm::Handle< std::vector< CaloParticle >> caloParticleHandle, std::vector< CaloParticle > const &cP, std::vector< size_t > const &cPIndices, std::vector< size_t > const &cPSelectedIndices, std::unordered_map< DetId, const HGCRecHit *> const &, std::map< double, double > cummatbudg, unsigned int layers, std::vector< int > thicknesses, const hgcal::RecoToSimCollection &recSimColl, const hgcal::SimToRecoCollection &simRecColl) const
Definition: HGVHistoProducerAlgo.cc:2176
HGVHistoProducerAlgo::getEta
double getEta(double eta) const
Definition: HGVHistoProducerAlgo.cc:3587
HCALHighEnergyHPDFilter_cfi.energy
energy
Definition: HCALHighEnergyHPDFilter_cfi.py:5
electrons_cff.bool
bool
Definition: electrons_cff.py:319
genParticles_cff.map
map
Definition: genParticles_cff.py:11
ticl::Trackster::vertices
std::vector< unsigned int > & vertices()
Definition: Trackster.h:57
DetId
Definition: DetId.h:17
mps_fire.end
end
Definition: mps_fire.py:242
edm::RefVector::size
size_type size() const
Size of the RefVector.
Definition: RefVector.h:102
HGVHistoProducerAlgo::fill_trackster_histos
void fill_trackster_histos(const Histograms &histograms, const int count, const ticl::TracksterCollection &Tracksters, const reco::CaloClusterCollection &layerClusters, const ticl::TracksterCollection &simTS, const ticl::TracksterCollection &simTS_fromCP, std::map< uint, std::vector< uint >> const &simTrackstersMap, std::vector< SimCluster > const &sC, const edm::ProductID &cPHandle_id, std::vector< CaloParticle > const &cP, std::vector< size_t > const &cPIndices, std::vector< size_t > const &cPSelectedIndices, std::unordered_map< DetId, const HGCRecHit *> const &, unsigned int layers) const
Definition: HGVHistoProducerAlgo.cc:3332
ticl::Trackster::seedID
const edm::ProductID & seedID() const
Definition: Trackster.h:125
HGVHistoProducerAlgo::recHitTools_
std::shared_ptr< hgcal::RecHitTools > recHitTools_
Definition: HGVHistoProducerAlgo.h:386
DetId::rawId
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:57
ticl::Trackster::vertex_multiplicity
std::vector< float > & vertex_multiplicity()
Definition: Trackster.h:58
HGVHistoProducerAlgo::tracksters_to_SimTracksters
void tracksters_to_SimTracksters(const Histograms &histograms, const int count, const ticl::TracksterCollection &Tracksters, const reco::CaloClusterCollection &layerClusters, const ticl::TracksterCollection &simTS, const validationType valType, const ticl::TracksterCollection &simTS_fromCP, std::map< uint, std::vector< uint >> const &simTrackstersMap, std::vector< SimCluster > const &sC, const edm::ProductID &cPHandle_id, std::vector< CaloParticle > const &cP, std::vector< size_t > const &cPIndices, std::vector< size_t > const &cPSelectedIndices, std::unordered_map< DetId, const HGCRecHit *> const &, unsigned int layers) const
Definition: HGVHistoProducerAlgo.cc:2516
dqm::implementation::IBooker::book2D
MonitorElement * book2D(TString const &name, TString const &title, int nchX, double lowX, double highX, int nchY, double lowY, double highY, FUNC onbooking=NOOP())
Definition: DQMStore.h:212
dqm::implementation::IBooker::bookInt
MonitorElement * bookInt(TString const &name, FUNC onbooking=NOOP())
Definition: DQMStore.h:73
ScoreCutSTStoTSPurDup_
const double ScoreCutSTStoTSPurDup_[]
Definition: HGVHistoProducerAlgo.cc:21
HGVHistoProducerAlgo::bookClusterHistos_ClusterLevel
void bookClusterHistos_ClusterLevel(DQMStore::IBooker &ibook, Histograms &histograms, unsigned int layers, std::vector< int > thicknesses, std::string pathtomatbudfile)
Definition: HGVHistoProducerAlgo.cc:600
testProducerWithPsetDescEmpty_cfi.y1
y1
Definition: testProducerWithPsetDescEmpty_cfi.py:29
ticl::Trackster::seedIndex
const int seedIndex() const
Definition: Trackster.h:126
position
static int position[264][3]
Definition: ReadPGInfo.cc:289
ticl::TracksterCollection
std::vector< Trackster > TracksterCollection
Definition: Trackster.h:203
reco::CaloCluster::eta
double eta() const
pseudorapidity of cluster centroid
Definition: CaloCluster.h:181
HGVHistoProducerAlgo::bookClusterHistos_CellLevel
void bookClusterHistos_CellLevel(DQMStore::IBooker &ibook, Histograms &histograms, unsigned int layers, std::vector< int > thicknesses)
Definition: HGVHistoProducerAlgo.cc:864
dqm::implementation::IBooker::book1D
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX, FUNC onbooking=NOOP())
Definition: DQMStore.h:98
HGVHistoProducerAlgo::layerClusters_to_CaloParticles
void layerClusters_to_CaloParticles(const Histograms &histograms, edm::Handle< reco::CaloClusterCollection > clusterHandle, const reco::CaloClusterCollection &clusters, edm::Handle< std::vector< CaloParticle >> caloParticleHandle, std::vector< CaloParticle > const &cP, std::vector< size_t > const &cPIndices, std::vector< size_t > const &cPSelectedIndices, std::unordered_map< DetId, const HGCRecHit *> const &, unsigned int layers, const hgcal::RecoToSimCollection &recSimColl, const hgcal::SimToRecoCollection &simRecColl) const
Definition: HGVHistoProducerAlgo.cc:1739
HGVHistoProducerAlgo::distance
double distance(const double x1, const double y1, const double x2, const double y2) const
Definition: HGVHistoProducerAlgo.cc:3555
h
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
testProducerWithPsetDescEmpty_cfi.x1
x1
Definition: testProducerWithPsetDescEmpty_cfi.py:33
hgcalTopologyTester_cfi.layers
layers
Definition: hgcalTopologyTester_cfi.py:8
Density
hgcal_clustering::Density Density
Definition: HGCalImagingAlgo.h:29
heppy_batch.val
val
Definition: heppy_batch.py:351
ScoreCutCPtoLC_
const double ScoreCutCPtoLC_
Definition: HGVHistoProducerAlgo.cc:17
CaloParticle::eta
float eta() const
Momentum pseudorapidity. Note this is taken from the simtrack before the calorimeter.
Definition: CaloParticle.h:142
FastTrackerRecHitCombiner_cfi.simHits
simHits
Definition: FastTrackerRecHitCombiner_cfi.py:5
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
HGCalValidator_cfi.simVertices
simVertices
Definition: HGCalValidator_cfi.py:66
eostools.move
def move(src, dest)
Definition: eostools.py:511
AlCaHLTBitMon_ParallelJobs.p
def p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
HGVHistoProducerAlgo::HGVHistoProducerAlgo
HGVHistoProducerAlgo(const edm::ParameterSet &pset)
Definition: HGVHistoProducerAlgo.cc:23
bsc_activity_cfg.clusters
clusters
Definition: bsc_activity_cfg.py:36
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:233
Generated for CMSSW Reference Manual by doxygen 1.8.14