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