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