CMS 3D CMS Logo

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