test
CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HcalRecHitsValidation.cc
Go to the documentation of this file.
4 
6  // DQM ROOT output
7  outputFile_ = conf.getUntrackedParameter<std::string>("outputFile", "myfile.root");
8 
9  if ( outputFile_.size() != 0 ) {
10  edm::LogInfo("OutputInfo") << " Hcal RecHit Task histograms will be saved to '" << outputFile_.c_str() << "'";
11  } else {
12  edm::LogInfo("OutputInfo") << " Hcal RecHit Task histograms will NOT be saved";
13  }
14 
15  nevtot = 0;
16 
17  hcalselector_ = conf.getUntrackedParameter<std::string>("hcalselector", "all");
18  ecalselector_ = conf.getUntrackedParameter<std::string>("ecalselector", "yes");
19  eventype_ = conf.getUntrackedParameter<std::string>("eventype", "single");
20  sign_ = conf.getUntrackedParameter<std::string>("sign", "*");
21  mc_ = conf.getUntrackedParameter<std::string>("mc", "yes");
22  famos_ = conf.getUntrackedParameter<bool>("Famos", false);
23  useAllHistos_ = conf.getUntrackedParameter<bool>("useAllHistos", false);
24 
25  //Collections
26  tok_hbhe_ = consumes<HBHERecHitCollection>(conf.getUntrackedParameter<edm::InputTag>("HBHERecHitCollectionLabel"));
27  tok_hf_ = consumes<HFRecHitCollection>(conf.getUntrackedParameter<edm::InputTag>("HFRecHitCollectionLabel"));
28  tok_ho_ = consumes<HORecHitCollection>(conf.getUntrackedParameter<edm::InputTag>("HORecHitCollectionLabel"));
29 
30  // register for data access
31  tok_evt_ = consumes<edm::HepMCProduct>(edm::InputTag("generatorSmeared"));
32  tok_EB_ = consumes<EBRecHitCollection>(edm::InputTag("ecalRecHit","EcalRecHitsEB"));
33  tok_EE_ = consumes<EERecHitCollection>(edm::InputTag("ecalRecHit","EcalRecHitsEE"));
34  tok_hh_ = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits","HcalHits"));
35 
36  // std::cout << "*** famos_ = " << famos_ << std::endl;
37 
38  subdet_ = 5;
39  if (hcalselector_ == "noise") subdet_ = 0;
40  if (hcalselector_ == "HB" ) subdet_ = 1;
41  if (hcalselector_ == "HE" ) subdet_ = 2;
42  if (hcalselector_ == "HO" ) subdet_ = 3;
43  if (hcalselector_ == "HF" ) subdet_ = 4;
44  if (hcalselector_ == "all" ) subdet_ = 5;
45  if (hcalselector_ == "ZS" ) subdet_ = 6;
46 
47  etype_ = 1;
48  if (eventype_ == "multi") etype_ = 2;
49 
50  iz = 1;
51  if(sign_ == "-") iz = -1;
52  if(sign_ == "*") iz = 0;
53 
54  imc = 1;
55  if(mc_ == "no") imc = 0;
56 
57 }
58 
59 
61 
63 {
64 
65  Char_t histo[200];
66 
67  ib.setCurrentFolder("HcalRecHitsV/HcalRecHitTask");
68 
69 
70 
71  // General counters (drawn)
72  sprintf (histo, "N_HB" );
73  Nhb = ib.book1D(histo, histo, 2600,0.,2600.);
74  sprintf (histo, "N_HE" );
75  Nhe = ib.book1D(histo, histo, 2600,0.,2600.);
76  sprintf (histo, "N_HO" );
77  Nho = ib.book1D(histo, histo, 2200,0.,2200.);
78  sprintf (histo, "N_HF" );
79  Nhf = ib.book1D(histo, histo, 1800,0., 1800.);
80 
81  // ZS
82  if(subdet_ == 6) {
83 
84  for (unsigned int i1 = 0; i1 < 82; i1++) {
85  for (unsigned int i2 = 0; i2 < 72; i2++) {
86  for (unsigned int i3 = 0; i3 < 4; i3++) {
87  for (unsigned int i4 = 0; i4 < 4; i4++) {
88  emap_min [i1][i2][i3][i4] = 99999.;
89  }
90  }
91  }
92  }
93 
94  //None of the ZS histos are drawn
95  if (useAllHistos_){
96  sprintf (histo, "ZSmin_map_depth1" );
97  map_depth1 = ib.book2D(histo, histo, 82, -41., 41., 72, 0., 72.);
98  sprintf (histo, "ZSmin_map_depth2" );
99  map_depth2 = ib.book2D(histo, histo, 82, -41., 41., 72, 0., 72.);
100  sprintf (histo, "ZSmin_map_depth3" );
101  map_depth3 = ib.book2D(histo, histo, 82, -41., 41., 72, 0., 72.);
102  sprintf (histo, "ZSmin_map_depth4" );
103  map_depth4 = ib.book2D(histo, histo, 82, -41., 41., 72, 0., 72.);
104 
105 
106  sprintf (histo, "ZS_Nreco_HB1" );
107  ZS_nHB1 = ib.book1D(histo, histo, 2500, 0., 2500.);
108  sprintf (histo, "ZS_Nreco_HB2" );
109  ZS_nHB2 = ib.book1D(histo, histo, 500, 0., 500.);
110  sprintf (histo, "ZS_Nreco_HE1" );
111  ZS_nHE1 = ib.book1D(histo, histo, 2000, 0., 2000.);
112  sprintf (histo, "ZS_Nreco_HE2" );
113  ZS_nHE2 = ib.book1D(histo, histo, 2000, 0., 2000.);
114  sprintf (histo, "ZS_Nreco_HE3" );
115  ZS_nHE3 = ib.book1D(histo, histo, 500, 0., 500.);
116  sprintf (histo, "ZS_Nreco_HO" );
117  ZS_nHO = ib.book1D(histo, histo, 2500, 0., 2500.);
118  sprintf (histo, "ZS_Nreco_HF1" );
119  ZS_nHF1 = ib.book1D(histo, histo, 1000, 0., 1000.);
120  sprintf (histo, "ZS_Nreco_HF2" );
121  ZS_nHF2 = ib.book1D(histo, histo, 1000, 0., 1000.);
122 
123  sprintf (histo, "ZSmin_simple1D_HB1" );
124  ZS_HB1 = ib.book1D(histo, histo,120, -2., 10.);
125  sprintf (histo, "ZSmin_simple1D_HB2" );
126  ZS_HB2 = ib.book1D(histo, histo,120, -2., 10.);
127  sprintf (histo, "ZSmin_simple1D_HE1" );
128  ZS_HE1 = ib.book1D(histo, histo,120, -2., 10.);
129  sprintf (histo, "ZSmin_simple1D_HE2" );
130  ZS_HE2 = ib.book1D(histo, histo,120, -2., 10.);
131  sprintf (histo, "ZSmin_simple1D_HE3" );
132  ZS_HE3 = ib.book1D(histo, histo,120, -2., 10.);
133  sprintf (histo, "ZSmin_simple1D_HO" );
134  ZS_HO = ib.book1D(histo, histo,120, -2., 10.);
135  sprintf (histo, "ZSmin_simple1D_HF1" );
136  ZS_HF1 = ib.book1D(histo, histo,200, -10., 10.);
137  sprintf (histo, "ZSmin_simple1D_HF2" );
138  ZS_HF2 = ib.book1D(histo, histo,200, -10., 10.);
139 
140  sprintf (histo, "ZSmin_sequential1D_HB1" );
141  ZS_seqHB1 = ib.book1D(histo, histo,2400, -1200., 1200.);
142  sprintf (histo, "ZSmin_sequential1D_HB2" );
143  ZS_seqHB2 = ib.book1D(histo, histo,2400, -1200., 1200.);
144  sprintf (histo, "ZSmin_sequential1D_HE1" );
145  ZS_seqHE1 = ib.book1D(histo, histo,4400, -2200., 2200.);
146  sprintf (histo, "ZSmin_sequential1D_HE2" );
147  ZS_seqHE2 = ib.book1D(histo, histo,4400, -2200., 2200.);
148  sprintf (histo, "ZSmin_sequential1D_HE3" );
149  ZS_seqHE3 = ib.book1D(histo, histo,4400, -2200., 2200.);
150  sprintf (histo, "ZSmin_sequential1D_HO" );
151  ZS_seqHO = ib.book1D(histo, histo,2400, -1200., 1200.);
152  sprintf (histo, "ZSmin_sequential1D_HF1" );
153  ZS_seqHF1 = ib.book1D(histo, histo,6000, -3000., 3000.);
154  sprintf (histo, "ZSmin_sequential1D_HF2" );
155  ZS_seqHF2 = ib.book1D(histo, histo,6000, -3000., 3000.);
156  }
157  }
158 
159  // ALL others, except ZS
160  else {
161 
162  sprintf (histo, "emap_depth1" );
163  emap_depth1 = ib.book2D(histo, histo, 84, -42., 42., 72, 0., 72.);
164  sprintf (histo, "emap_depth2" );
165  emap_depth2 = ib.book2D(histo, histo, 84, -42., 42., 72, 0., 72.);
166  sprintf (histo, "emap_depth3" );
167  emap_depth3 = ib.book2D(histo, histo, 84, -42., 42., 72, 0., 72.);
168  sprintf (histo, "emap_depth4" );
169  emap_depth4 = ib.book2D(histo, histo, 84, -42., 42., 72, 0., 72.);
170 
171  if (useAllHistos_){
172 
173  if (ecalselector_ == "yes") {
174  sprintf (histo, "map_ecal" );
175  map_ecal = ib.book2D(histo, histo, 70, -3.045, 3.045, 72, -3.1415926536, 3.1415926536);
176  }
177  }
178 
179  //The mean energy histos are drawn, but not the RMS or emean seq
180  sprintf (histo, "emean_vs_ieta_HB1" );
181  emean_vs_ieta_HB1 = ib.bookProfile(histo, histo, 82, -41., 41., 2010, -10., 2000., "s");
182  sprintf (histo, "emean_vs_ieta_HB2" );
183  emean_vs_ieta_HB2 = ib.bookProfile(histo, histo, 82, -41., 41., 2010, -10., 2000., "s");
184  sprintf (histo, "emean_vs_ieta_HE1" );
185  emean_vs_ieta_HE1 = ib.bookProfile(histo, histo, 82, -41., 41., 2010, -10. ,2000., "s" );
186  sprintf (histo, "emean_vs_ieta_HE2" );
187  emean_vs_ieta_HE2 = ib.bookProfile(histo, histo, 82, -41., 41., 2010, -10., 2000., "s");
188  sprintf (histo, "emean_vs_ieta_HE3" );
189  emean_vs_ieta_HE3 = ib.bookProfile(histo, histo, 82, -41., 41., 2010, -10., 2000., "s" );
190  sprintf (histo, "emean_vs_ieta_HO" );
191  emean_vs_ieta_HO = ib.bookProfile(histo, histo, 82, -41., 41., 2010, -10., 2000., "s" );
192  sprintf (histo, "emean_vs_ieta_HF1" );
193  emean_vs_ieta_HF1 = ib.bookProfile(histo, histo, 82, -41., 41., 2010, -10., 2000., "s" );
194  sprintf (histo, "emean_vs_ieta_HF2" );
195  emean_vs_ieta_HF2 = ib.bookProfile(histo, histo, 82, -41., 41., 2010, -10., 2000., "s" );
196 
197  if (useAllHistos_){
198  sprintf (histo, "RMS_vs_ieta_HB1" );
199  RMS_vs_ieta_HB1 = ib.book1D(histo, histo, 82, -41., 41.);
200  sprintf (histo, "RMS_vs_ieta_HB2" );
201  RMS_vs_ieta_HB2 = ib.book1D(histo, histo, 82, -41., 41.);
202  sprintf (histo, "RMS_vs_ieta_HE1" );
203  RMS_vs_ieta_HE1 = ib.book1D(histo, histo, 82, -41., 41.);
204  sprintf (histo, "RMS_vs_ieta_HE2" );
205  RMS_vs_ieta_HE2 = ib.book1D(histo, histo, 82, -41., 41.);
206  sprintf (histo, "RMS_vs_ieta_HE3" );
207  RMS_vs_ieta_HE3 = ib.book1D(histo, histo, 82, -41., 41.);
208  sprintf (histo, "RMS_vs_ieta_HO" );
209  RMS_vs_ieta_HO = ib.book1D(histo, histo, 82, -41., 41.);
210  sprintf (histo, "RMS_vs_ieta_HF1" );
211  RMS_vs_ieta_HF1 = ib.book1D(histo, histo, 82, -41., 41.);
212  sprintf (histo, "RMS_vs_ieta_HF2" );
213  RMS_vs_ieta_HF2 = ib.book1D(histo, histo, 82, -41., 41.);
214 
215  // Sequential emean and RMS
216  sprintf (histo, "emean_seq_HB1" );
217  emean_seqHB1 = ib.bookProfile(histo, histo, 2400, -1200., 1200., 2010, -10., 2000., "s" );
218  sprintf (histo, "emean_seq_HB2" );
219  emean_seqHB2 = ib.bookProfile(histo, histo, 2400, -1200., 1200., 2010, -10., 2000., "s" );
220  sprintf (histo, "emean_seq_HE1" );
221  emean_seqHE1 = ib.bookProfile(histo, histo, 4400, -2200., 2200., 2010, -10., 2000., "s" );
222  sprintf (histo, "emean_seq_HE2" );
223  emean_seqHE2 = ib.bookProfile(histo, histo, 4400, -2200., 2200., 2010, -10., 2000., "s" );
224  sprintf (histo, "emean_seq_HE3" );
225  emean_seqHE3 = ib.bookProfile(histo, histo, 4400, -2200., 2200., 2010, -10., 2000., "s" );
226  sprintf (histo, "emean_seq_HO" );
227  emean_seqHO = ib.bookProfile(histo, histo, 2400, -1200., 1200., 2010, -10., 2000., "s" );
228  sprintf (histo, "emean_seq_HF1" );
229  emean_seqHF1 = ib.bookProfile(histo, histo, 6000, -3000., 3000., 2010, -10., 2000., "s" );
230  sprintf (histo, "emean_seq_HF2" );
231  emean_seqHF2 = ib.bookProfile(histo, histo, 6000, -3000., 3000., 2010, -10., 2000., "s" );
232 
233  sprintf (histo, "RMS_seq_HB1" );
234  RMS_seq_HB1 = ib.book1D(histo, histo, 2400, -1200., 1200.);
235  sprintf (histo, "RMS_seq_HB2" );
236  RMS_seq_HB2 = ib.book1D(histo, histo, 2400, -1200., 1200.);
237  sprintf (histo, "RMS_seq_HE1" );
238  RMS_seq_HE1 = ib.book1D(histo, histo, 4400, -2200., 2200.);
239  sprintf (histo, "RMS_seq_HE2" );
240  RMS_seq_HE2 = ib.book1D(histo, histo, 4400, -2200., 2200.);
241  sprintf (histo, "RMS_seq_HE3" );
242  RMS_seq_HE3 = ib.book1D(histo, histo, 4400, -2200., 2200.);
243  sprintf (histo, "RMS_seq_HO" );
244  RMS_seq_HO = ib.book1D(histo, histo, 2400, -1200., 1200.);
245  sprintf (histo, "RMS_seq_HF1" );
246  RMS_seq_HF1 = ib.book1D(histo, histo, 6000, -3000., 3000.);
247  sprintf (histo, "RMS_seq_HF2" );
248  RMS_seq_HF2 = ib.book1D(histo, histo, 6000, -3000., 3000.);
249  }
250  // Occupancy
251  //The only occupancy histos drawn are occupancy vs. ieta
252  //but the maps are needed because this is where the latter are filled from
253  sprintf (histo, "occupancy_map_HB1" );
254  occupancy_map_HB1 = ib.book2D(histo, histo, 82, -41., 41., 72, 0., 72.);
255  sprintf (histo, "occupancy_map_HB2" );
256  occupancy_map_HB2 = ib.book2D(histo, histo, 82, -41., 41., 72, 0., 72.);
257  sprintf (histo, "occupancy_map_HE1" );
258  occupancy_map_HE1 = ib.book2D(histo, histo, 82, -41., 41., 72, 0., 72.);
259  sprintf (histo, "occupancy_map_HE2" );
260  occupancy_map_HE2 = ib.book2D(histo, histo, 82, -41., 41., 72, 0., 72.);
261  sprintf (histo, "occupancy_map_HE3" );
262  occupancy_map_HE3 = ib.book2D(histo, histo, 82, -41., 41., 72, 0., 72.);
263  sprintf (histo, "occupancy_map_HO" );
264  occupancy_map_HO = ib.book2D(histo, histo, 82, -41., 41., 72, 0., 72.);
265  sprintf (histo, "occupancy_map_HF1" );
266  occupancy_map_HF1 = ib.book2D(histo, histo, 82, -41., 41., 72, 0., 72.);
267  sprintf (histo, "occupancy_map_HF2" );
268  occupancy_map_HF2 = ib.book2D(histo, histo, 82, -41., 41., 72, 0., 72.);
269 
270  //These are drawn
271  sprintf (histo, "occupancy_vs_ieta_HB1" );
272  occupancy_vs_ieta_HB1 = ib.book1D(histo, histo, 82, -41., 41.);
273  sprintf (histo, "occupancy_vs_ieta_HB2" );
274  occupancy_vs_ieta_HB2 = ib.book1D(histo, histo, 82, -41., 41.);
275  sprintf (histo, "occupancy_vs_ieta_HE1" );
276  occupancy_vs_ieta_HE1 = ib.book1D(histo, histo, 82, -41., 41.);
277  sprintf (histo, "occupancy_vs_ieta_HE2" );
278  occupancy_vs_ieta_HE2 = ib.book1D(histo, histo, 82, -41., 41.);
279  sprintf (histo, "occupancy_vs_ieta_HE3" );
280  occupancy_vs_ieta_HE3 = ib.book1D(histo, histo, 82, -41., 41.);
281  sprintf (histo, "occupancy_vs_ieta_HO" );
282  occupancy_vs_ieta_HO = ib.book1D(histo, histo, 82, -41., 41.);
283  sprintf (histo, "occupancy_vs_ieta_HF1" );
284  occupancy_vs_ieta_HF1 = ib.book1D(histo, histo, 82, -41., 41.);
285  sprintf (histo, "occupancy_vs_ieta_HF2" );
286  occupancy_vs_ieta_HF2 = ib.book1D(histo, histo, 82, -41., 41.);
287 
288  //These are not
289  if (useAllHistos_){
290  sprintf (histo, "occ_sequential1D_HB1" );
291  occupancy_seqHB1 = ib.book1D(histo, histo,2400, -1200., 1200.);
292  sprintf (histo, "occ_sequential1D_HB2" );
293  occupancy_seqHB2 = ib.book1D(histo, histo,2400, -1200., 1200.);
294  sprintf (histo, "occ_sequential1D_HE1" );
295  occupancy_seqHE1 = ib.book1D(histo, histo,4400, -2200., 2200.);
296  sprintf (histo, "occ_sequential1D_HE2" );
297  occupancy_seqHE2 = ib.book1D(histo, histo,4400, -2200., 2200.);
298  sprintf (histo, "occ_sequential1D_HE3" );
299  occupancy_seqHE3 = ib.book1D(histo, histo,4400, -2200., 2200.);
300  sprintf (histo, "occ_sequential1D_HO" );
301  occupancy_seqHO = ib.book1D(histo, histo,2400, -1200., 1200.);
302  sprintf (histo, "occ_sequential1D_HF1" );
303  occupancy_seqHF1 = ib.book1D(histo, histo,6000, -3000., 3000.);
304  sprintf (histo, "occ_sequential1D_HF2" );
305  occupancy_seqHF2 = ib.book1D(histo, histo,6000, -3000., 3000.);
306  }
307 
308  //All status word histos except HF67 are drawn
309  sprintf (histo, "HcalRecHitTask_RecHit_StatusWord_HB" ) ;
310  RecHit_StatusWord_HB = ib.book1D(histo, histo, 32 , -0.5, 31.5);
311 
312  sprintf (histo, "HcalRecHitTask_RecHit_StatusWord_HE" ) ;
313  RecHit_StatusWord_HE = ib.book1D(histo, histo, 32 , -0.5, 31.5);
314 
315  sprintf (histo, "HcalRecHitTask_RecHit_StatusWord_HF" ) ;
316  RecHit_StatusWord_HF = ib.book1D(histo, histo, 32 , -0.5, 31.5);
317 
318  if (useAllHistos_){
319  sprintf (histo, "HcalRecHitTask_RecHit_StatusWord_HF67" ) ;
320  RecHit_StatusWord_HF67 = ib.book1D(histo, histo, 3 , 0.5, 3.5);
321  }
322  sprintf (histo, "HcalRecHitTask_RecHit_StatusWord_HO" ) ;
323  RecHit_StatusWord_HO = ib.book1D(histo, histo, 32 , -0.5, 31.5);
324 
325  //Aux status word histos
326  sprintf (histo, "HcalRecHitTask_RecHit_Aux_StatusWord_HB" ) ;
327  RecHit_Aux_StatusWord_HB = ib.book1D(histo, histo, 32 , -0.5, 31.5);
328 
329  sprintf (histo, "HcalRecHitTask_RecHit_Aux_StatusWord_HE" ) ;
330  RecHit_Aux_StatusWord_HE = ib.book1D(histo, histo, 32 , -0.5, 31.5);
331 
332  sprintf (histo, "HcalRecHitTask_RecHit_Aux_StatusWord_HF" ) ;
333  RecHit_Aux_StatusWord_HF = ib.book1D(histo, histo, 32 , -0.5, 31.5);
334 
335  sprintf (histo, "HcalRecHitTask_RecHit_Aux_StatusWord_HO" ) ;
336  RecHit_Aux_StatusWord_HO = ib.book1D(histo, histo, 32 , -0.5, 31.5);
337 
338  //Status word correlation plots
339  sprintf (histo, "HcalRecHitTask_RecHit_StatusWordCorr_HB");
340  RecHit_StatusWordCorr_HB = ib.book2D(histo, histo, 2, -0.5, 1.5, 2, -0.5, 1.5);
341 
342  sprintf (histo, "HcalRecHitTask_RecHit_StatusWordCorr_HE");
343  RecHit_StatusWordCorr_HE = ib.book2D(histo, histo, 2, -0.5, 1.5, 2, -0.5, 1.5);
344  //These are not drawn
345  if(imc !=0 && useAllHistos_) {
346  sprintf (histo, "map_econe_depth1" );
348  ib.book2D(histo, histo, 520, -5.2, 5.2, 72, -3.1415926536, 3.1415926536);
349  sprintf (histo, "map_econe_depth2" );
351  ib.book2D(histo, histo, 520, -5.2, 5.2, 72, -3.1415926536, 3.1415926536);
352  sprintf (histo, "map_econe_depth3" );
354  ib.book2D(histo, histo, 520, -5.2, 5.2, 72, -3.1415926536, 3.1415926536);
355  sprintf (histo, "map_econe_depth4" );
357  ib.book2D(histo, histo, 520, -5.2, 5.2, 72, -3.1415926536, 3.1415926536);
358  }
359  } // end-of (subdet_ =! 6)
360 
361  //======================= Now various cases one by one ===================
362 
363  //Histograms drawn for single pion scan
364  if(subdet_ != 0 && imc != 0) { // just not for noise
365  sprintf (histo, "HcalRecHitTask_En_rechits_cone_profile_vs_ieta_all_depths");
366  meEnConeEtaProfile = ib.bookProfile(histo, histo, 82, -41., 41., 2100, -100., 2000.);
367 
368  sprintf (histo, "HcalRecHitTask_En_rechits_cone_profile_vs_ieta_all_depths_E");
369  meEnConeEtaProfile_E = ib.bookProfile(histo, histo, 82, -41., 41., 2100, -100., 2000.);
370 
371  sprintf (histo, "HcalRecHitTask_En_rechits_cone_profile_vs_ieta_all_depths_EH");
372  meEnConeEtaProfile_EH = ib.bookProfile(histo, histo, 82, -41., 41., 2100, -100., 2000.);
373  }
374  //The other cone profile, delta ieta/phi and noise histos are not drawn
375  if (useAllHistos_){
376  if(subdet_ != 0 && imc != 0) { // just not for noise
377 
378  // meEnConeEtaProfiel_depth1->Fill(eta_RecHit, HcalCone_d1);
379 
380  sprintf (histo, "HcalRecHitTask_En_rechits_cone_profile_vs_ieta_depth1");
381  meEnConeEtaProfile_depth1 = ib.bookProfile(histo, histo, 82, -41., 41., 2100, -100., 2000.);
382 
383  sprintf (histo, "HcalRecHitTask_En_rechits_cone_profile_vs_ieta_depth2");
384  meEnConeEtaProfile_depth2 = ib.bookProfile(histo, histo, 82, -41., 41., 2100, -100., 2000.);
385 
386  sprintf (histo, "HcalRecHitTask_En_rechits_cone_profile_vs_ieta_depth3");
387  meEnConeEtaProfile_depth3 = ib.bookProfile(histo, histo, 82, -41., 41., 2100, -100., 2000.);
388 
389  sprintf (histo, "HcalRecHitTask_En_rechits_cone_profile_vs_ieta_depth4");
390  meEnConeEtaProfile_depth4 = ib.bookProfile(histo, histo, 82, -41., 41., 2100, -100., 2000.);
391 
392  }
393 
394  if(etype_ == 1 && subdet_ != 0) { // single part., not for noise
395 
396  sprintf (histo, "Delta_phi_cluster-MC");
397  meDeltaPhi = ib.book2D(histo, histo, 520, -5.2, 5.2, 60, -0.6, 0.6);
398 
399  sprintf (histo, "Delta_eta_cluster-MC");
400  meDeltaEta = ib.book2D(histo, histo, 520, -5.2, 5.2, 60, -0.6, 0.6);
401 
402  sprintf (histo, "Delta_phi_simcluster-MC");
403  meDeltaPhiS = ib.book2D(histo, histo, 520, -5.2, 5.2, 60, -0.6, 0.6);
404 
405  sprintf (histo, "Delta_eta_simcluster-MC");
406  meDeltaEtaS = ib.book2D(histo, histo, 520, -5.2, 5.2, 60, -0.6, 0.6);
407  }
408  // NOISE-specific
409 
410  if (hcalselector_ == "noise" ){
411 
412  sprintf (histo, "e_hb" ) ;
413  e_hb = ib.book1D(histo, histo,1000, -5., 5.);
414  sprintf (histo, "e_he" ) ;
415  e_he = ib.book1D(histo, histo,1000, -5., 5.);
416  sprintf (histo, "e_ho" ) ;
417  e_ho = ib.book1D(histo, histo,1000, -5., 5.);
418  sprintf (histo, "e_hfl" ) ;
419  e_hfl = ib.book1D(histo, histo,2000, -10., 10.);
420  sprintf (histo, "e_hfs" ) ;
421  e_hfs = ib.book1D(histo, histo,2000, -10., 10.);
422  }
423  }
424  // ************** HB **********************************
425  if (subdet_ == 1 || subdet_ == 5 ){
426 
427  //Only severity level, energy of rechits and overall HB timing histos are drawn
428  if (useAllHistos_){
429  if(etype_ == 1 && subdet_ == 1 ) {
430  if(imc != 0) {
431  sprintf (histo, "HcalRecHitTask_number_of_rechits_in_cone_HB" ) ;
432  meNumRecHitsConeHB = ib.book1D(histo, histo, 100, 0., 100.);
433 
434  sprintf (histo, "HcalRecHitTask_sum_of_rechits_energy_in_cone_HB" ) ;
435  meSumRecHitsEnergyConeHB = ib.book1D(histo,histo, 60 ,-20., 280.);
436  }
437 
438  sprintf (histo, "HcalRecHitTask_number_of_rechits_above_1GeV_HB");
439  meNumRecHitsThreshHB = ib.book1D(histo, histo, 30, 0., 30.);
440 
441  sprintf (histo, "HcalRecHitTask_sum_of_rechits_energy_HB" ) ;
442  meSumRecHitsEnergyHB = ib.book1D(histo,histo, 60 , -20., 280.);
443 
444  if (ecalselector_ == "yes") {
445  if(imc != 0) {
446  sprintf (histo, "HcalRecHitTask_number_of_ecalrechits_in_cone_HB");
447  meNumEcalRecHitsConeHB = ib.book1D(histo, histo, 300, 0., 300.);
448  sprintf (histo, "HcalRecHitTask_energy_ecal_plus_hcal_in_cone_HB");
449  meEcalHcalEnergyConeHB = ib.book1D(histo,histo, 60 , -20., 280.);
450  }
451 
452  sprintf (histo, "HcalRecHitTask_energy_hcal_vs_ecal_HB");
453  meEnergyHcalVsEcalHB = ib.book2D(histo, histo, 300, 0., 150., 300, 0., 150.);
454  sprintf (histo, "HcalRecHitTask_energy_ecal_plus_hcal_HB" ) ;
455  meEcalHcalEnergyHB = ib.book1D(histo,histo, 60 , -20., 280.);
456  }
457  }
458  }
459 
460  sprintf(histo, "HcalRecHitTask_severityLevel_HB");
461  sevLvl_HB = ib.book1D(histo, histo, 25, -0.5, 24.5);
462 
463  sprintf (histo, "HcalRecHitTask_energy_of_rechits_HB" ) ;
464  meRecHitsEnergyHB = ib.book1D(histo, histo, 2010 , -10. , 2000.);
465 
466  sprintf (histo, "HcalRecHitTask_timing_HB" ) ;
467  meTimeHB = ib.book1D(histo, histo, 70, -48., 92.);
468 
469  //High, medium and low histograms to reduce RAM usage
470  sprintf (histo, "HcalRecHitTask_timing_vs_energy_Low_HB" ) ;
471  meTE_Low_HB = ib.book2D(histo, histo, 50, -5., 45., 70, -48., 92.);
472 
473  sprintf (histo, "HcalRecHitTask_timing_vs_energy_HB" ) ;
474  meTE_HB = ib.book2D(histo, histo, 150, -5., 295., 70, -48., 92.);
475 
476  sprintf (histo, "HcalRecHitTask_timing_vs_energy_High_HB" ) ;
477  meTE_High_HB = ib.book2D(histo, histo, 150, -5., 2995., 70, -48., 92.);
478 
479  sprintf (histo, "HcalRecHitTask_timing_vs_energy_profile_Low_HB" ) ;
480  meTEprofileHB_Low = ib.bookProfile(histo, histo, 50, -5., 45., 70, -48., 92.);
481 
482  sprintf (histo, "HcalRecHitTask_timing_vs_energy_profile_HB" ) ;
483  meTEprofileHB = ib.bookProfile(histo, histo, 150, -5., 295., 70, -48., 92.);
484 
485  sprintf (histo, "HcalRecHitTask_timing_vs_energy_profile_High_HB" ) ;
486  meTEprofileHB_High = ib.bookProfile(histo, histo, 150, -5., 2995., 70, -48., 92.);
487 
488  //Timing by depth and rechits vs simhits are not drawn
489  if (useAllHistos_){
490  sprintf (histo, "HcalRecHitTask_timing_vs_energy_HB_depth1" ) ;
491  meTE_HB1 = ib.book2D(histo, histo, 3000, -5., 2995., 70, -48., 92.);
492 
493  sprintf (histo, "HcalRecHitTask_timing_vs_energy_HB_depth2" ) ;
494  meTE_HB2 = ib.book2D(histo, histo, 3000, -5., 2995., 70, -48., 92.);
495 
496  if(imc != 0) {
497  sprintf (histo, "HcalRecHitTask_energy_rechits_vs_simhits_HB");
498  meRecHitSimHitHB = ib.book2D(histo, histo, 120, 0., 1.2, 300, 0., 150.);
499  sprintf (histo, "HcalRecHitTask_energy_rechits_vs_simhits_profile_HB");
500  meRecHitSimHitProfileHB = ib.bookProfile(histo, histo, 120, 0., 1.2, 500, 0., 500.);
501  }
502  }
503  }
504 
505  // ********************** HE ************************************
506  if ( subdet_ == 2 || subdet_ == 5 ){
507 
508  //None of these are drawn
509  if (useAllHistos_){
510  if(etype_ == 1 && subdet_ == 2 ) {
511 
512  if(imc != 0) {
513  sprintf (histo, "HcalRecHitTask_number_of_rechits_in_cone_HE" ) ;
514  meNumRecHitsConeHE = ib.book1D(histo, histo, 100, 0., 100.);
515 
516  sprintf (histo, "HcalRecHitTask_sum_of_rechits_energy_in_cone_HE" ) ;
517  meSumRecHitsEnergyConeHE = ib.book1D(histo,histo, 60 ,-20., 280.);
518  }
519 
520  sprintf (histo, "HcalRecHitTask_number_of_rechits_above_1GeV_HE");
521  meNumRecHitsThreshHE = ib.book1D(histo, histo, 30, 0., 30.);
522 
523  sprintf (histo, "HcalRecHitTask_sum_of_rechits_energy_HE" ) ;
524  meSumRecHitsEnergyHE = ib.book1D(histo,histo, 60 , -20., 280.);
525 
526  if (ecalselector_ == "yes") {
527  sprintf (histo, "HcalRecHitTask_energy_ecal_plus_hcal_HE" ) ;
528  meEcalHcalEnergyHE = ib.book1D(histo,histo, 80, -20., 380.);
529 
530  sprintf (histo, "HcalRecHitTask_energy_hcal_vs_ecal_HE");
531  meEnergyHcalVsEcalHE = ib.book2D(histo, histo, 300, 0., 150., 300, 0., 150.);
532  if(imc != 0) {
533  sprintf (histo, "HcalRecHitTask_number_of_ecalrechits_in_cone_HE");
534  meNumEcalRecHitsConeHE = ib.book1D(histo, histo, 300, 0., 300.);
535  sprintf (histo, "HcalRecHitTask_energy_ecal_plus_hcal_in_cone_HE");
536  meEcalHcalEnergyConeHE = ib.book1D(histo,histo, 60,-20., 280.);
537  }
538  }
539  }
540  }
541 
542  //Only severity level, energy of rechits and overall HB timing histos are drawn
543  sprintf(histo, "HcalRecHitTask_severityLevel_HE");
544  sevLvl_HE = ib.book1D(histo, histo, 25, -0.5, 24.5);
545 
546  sprintf (histo, "HcalRecHitTask_energy_of_rechits_HE" ) ;
547  meRecHitsEnergyHE = ib.book1D(histo, histo, 2010, -10., 2000.);
548 
549  sprintf (histo, "HcalRecHitTask_timing_HE" ) ;
550  meTimeHE = ib.book1D(histo, histo, 70, -48., 92.);
551 
552  sprintf (histo, "HcalRecHitTask_timing_vs_energy_Low_HE" ) ;
553  meTE_Low_HE = ib.book2D(histo, histo, 80, -5., 75., 70, -48., 92.);
554 
555  sprintf (histo, "HcalRecHitTask_timing_vs_energy_HE" ) ;
556  meTE_HE = ib.book2D(histo, histo, 200, -5., 2995., 70, -48., 92.);
557 
558  sprintf (histo, "HcalRecHitTask_timing_vs_energy_profile_Low_HE" ) ;
559  meTEprofileHE_Low = ib.bookProfile(histo, histo, 80, -5., 75., 70, -48., 92.);
560 
561  sprintf (histo, "HcalRecHitTask_timing_vs_energy_profile_HE" ) ;
562  meTEprofileHE = ib.bookProfile(histo, histo, 200, -5., 2995., 70, -48., 92.);
563 
564  //Timing by depth and rechits vs simhits are not drawn
565  if (useAllHistos_){
566  sprintf (histo, "HcalRecHitTask_timing_vs_energy_HE_depth1" ) ;
567  meTE_HE1 = ib.book2D(histo, histo, 1000, -5., 995., 70, -48., 92.);
568 
569  sprintf (histo, "HcalRecHitTask_timing_vs_energy_HE_depth2" ) ;
570  meTE_HE2 = ib.book2D(histo, histo, 1000, -5., 995., 70, -48., 92.);
571 
572  if(imc != 0) {
573  sprintf (histo, "HcalRecHitTask_energy_rechits_vs_simhits_HE");
574  meRecHitSimHitHE = ib.book2D(histo, histo, 120, 0., 0.6, 300, 0., 150.);
575  sprintf (histo, "HcalRecHitTask_energy_rechits_vs_simhits_profile_HE");
576  meRecHitSimHitProfileHE = ib.bookProfile(histo, histo, 120, 0., 0.6, 500, 0., 500.);
577  }
578  }
579 
580  }
581 
582  // ************** HO ****************************************
583  if ( subdet_ == 3 || subdet_ == 5 ){
584 
585  //Only severity level, energy of rechits and overall HB timing histos are drawn
586  if (useAllHistos_){
587  if(etype_ == 1 && subdet_ == 3) {
588  if (imc != 0) {
589  sprintf (histo, "HcalRecHitTask_number_of_rechits_in_cone_HO" ) ;
590  meNumRecHitsConeHO = ib.book1D(histo, histo, 100, 0 , 100.);
591 
592  sprintf (histo, "HcalRecHitTask_sum_of_rechits_energy_in_cone_HO" ) ;
593  meSumRecHitsEnergyConeHO = ib.book1D(histo,histo, 80 ,-20., 380.);
594  }
595 
596  sprintf (histo, "HcalRecHitTask_number_of_rechits_above_1GeV_HO");
597  meNumRecHitsThreshHO = ib.book1D(histo, histo, 100, 0., 100.);
598 
599  sprintf (histo, "HcalRecHitTask_sum_of_rechits_energy_HO" ) ;
600  meSumRecHitsEnergyHO = ib.book1D(histo,histo, 80 , -20., 380.);
601  }
602  }
603 
604  sprintf(histo, "HcalRecHitTask_severityLevel_HO");
605  sevLvl_HO = ib.book1D(histo, histo, 25, -0.5, 24.5);
606 
607  sprintf (histo, "HcalRecHitTask_energy_of_rechits_HO" ) ;
608  meRecHitsEnergyHO = ib.book1D(histo, histo, 2010 , -10. , 2000.);
609 
610  sprintf (histo, "HcalRecHitTask_timing_HO" ) ;
611  meTimeHO = ib.book1D(histo, histo, 70, -48., 92.);
612 
613  sprintf (histo, "HcalRecHitTask_timing_vs_energy_HO" ) ;
614  meTE_HO= ib.book2D(histo, histo, 60, -5., 55., 70, -48., 92.);
615 
616  sprintf (histo, "HcalRecHitTask_timing_vs_energy_High_HO" ) ;
617  meTE_High_HO= ib.book2D(histo, histo, 100, -5., 995., 70, -48., 92.);
618 
619  sprintf (histo, "HcalRecHitTask_timing_vs_energy_profile_HO" ) ;
620  meTEprofileHO = ib.bookProfile(histo, histo, 60, -5., 55., 70, -48., 92.);
621 
622  sprintf (histo, "HcalRecHitTask_timing_vs_energy_profile_High_HO" ) ;
623  meTEprofileHO_High = ib.bookProfile(histo, histo, 100, -5., 995., 70, -48., 92.);
624 
625  //Rechits vs simhits are not drawn
626  if (useAllHistos_){
627  if(imc != 0) {
628  sprintf (histo, "HcalRecHitTask_energy_rechits_vs_simhits_HO");
629  meRecHitSimHitHO = ib.book2D(histo, histo, 150, 0., 1.5, 350, 0., 350.);
630  sprintf (histo, "HcalRecHitTask_energy_rechits_vs_simhits_profile_HO");
631  meRecHitSimHitProfileHO = ib.bookProfile(histo, histo, 150, 0., 1.5, 500, 0., 500.);
632  }
633  }
634  }
635 
636  // ********************** HF ************************************
637  if ( subdet_ == 4 || subdet_ == 5 ){
638 
639  //Only severity level, energy of rechits and overall HB timing histos are drawn
640  if (useAllHistos_){
641  if(etype_ == 1 && subdet_ == 4) {
642 
643  if(imc != 0) {
644  sprintf (histo, "HcalRecHitTask_number_of_rechits_in_cone_HF" ) ;
645  meNumRecHitsConeHF = ib.book1D(histo, histo, 30, 0 , 30.);
646 
647  sprintf (histo, "HcalRecHitTask_sum_of_rechits_energy_in_cone_HF" ) ;
648  meSumRecHitsEnergyConeHF = ib.book1D(histo,histo,100, -20., 180.);
649 
650  sprintf (histo, "HcalRecHitTask_sum_of_rechits_energy_in_cone_HFL" );
651  meSumRecHitsEnergyConeHFL = ib.book1D(histo,histo,100,-20., 180.);
652 
653  sprintf (histo, "HcalRecHitTask_sum_of_rechits_energy_in_cone_HFS");
654  meSumRecHitsEnergyConeHFS = ib.book1D(histo,histo,100,-20., 180.);
655  }
656  sprintf (histo, "HcalRecHitTask_sum_of_rechits_energy_HF" ) ;
657  meSumRecHitsEnergyHF = ib.book1D(histo,histo, 80 , -20., 380.);
658  }
659  }
660 
661  sprintf(histo, "HcalRecHitTask_severityLevel_HF");
662  sevLvl_HF = ib.book1D(histo, histo, 25, -0.5, 24.5);
663 
664  sprintf (histo, "HcalRecHitTask_energy_of_rechits_HF" ) ;
665  meRecHitsEnergyHF = ib.book1D(histo, histo, 2010 , -10. , 2000.);
666 
667  sprintf (histo, "HcalRecHitTask_timing_HF" ) ;
668  meTimeHF = ib.book1D(histo, histo, 70, -48., 92.);
669 
670  sprintf (histo, "HcalRecHitTask_timing_vs_energy_Low_HF" ) ;
671  meTE_Low_HF = ib.book2D(histo, histo, 100, -5., 195., 70, -48., 92.);
672 
673  sprintf (histo, "HcalRecHitTask_timing_vs_energy_HF" ) ;
674  meTE_HF = ib.book2D(histo, histo, 200, -5., 995., 70, -48., 92.);
675 
676  sprintf (histo, "HcalRecHitTask_timing_vs_energy_profile_Low_HF" ) ;
677  meTEprofileHF_Low = ib.bookProfile(histo, histo, 100, -5., 195., 70, -48., 92.);
678 
679  sprintf (histo, "HcalRecHitTask_timing_vs_energy_profile_HF" ) ;
680  meTEprofileHF = ib.bookProfile(histo, histo, 200, -5., 995., 70, -48., 92.);
681 
682  //Timing by L/S and rechits vs simhits are not drawn
683  if (useAllHistos_){
684  sprintf (histo, "HcalRecHitTask_timing_vs_energy_HFL" ) ;
685  meTE_HFL = ib.book2D(histo, histo, 1000, -5., 995., 70, -48., 92.);
686 
687  sprintf (histo, "HcalRecHitTask_timing_vs_energy_HFS" ) ;
688  meTE_HFS = ib.book2D(histo, histo, 1000, -5., 995., 70, -48., 92.);
689 
690  if(imc != 0) {
691  sprintf (histo, "HcalRecHitTask_energy_rechits_vs_simhits_HF");
692  meRecHitSimHitHF = ib.book2D(histo, histo, 50, 0., 50., 150, 0., 150.);
693  sprintf (histo, "HcalRecHitTask_energy_rechits_vs_simhits_HFL");
694  meRecHitSimHitHFL = ib.book2D(histo, histo, 50, 0., 50., 150, 0., 150.);
695  sprintf (histo, "HcalRecHitTask_energy_rechits_vs_simhits_HFS");
696  meRecHitSimHitHFS = ib.book2D(histo, histo, 50, 0., 50., 150, 0., 150.);
697  sprintf (histo, "HcalRecHitTask_energy_rechits_vs_simhits_profile_HF");
698  meRecHitSimHitProfileHF = ib.bookProfile(histo, histo, 50, 0., 50., 500, 0., 500.);
699  sprintf (histo, "HcalRecHitTask_energy_rechits_vs_simhits_profile_HFL");
700  meRecHitSimHitProfileHFL = ib.bookProfile(histo, histo, 50, 0., 50., 500, 0., 500.);
701  sprintf (histo, "HcalRecHitTask_energy_rechits_vs_simhits_profile_HFS");
702  meRecHitSimHitProfileHFS = ib.bookProfile(histo, histo, 50, 0., 50., 500, 0., 500.);
703  }
704  }
705  }
706 
707 }
708 
710 
711  using namespace edm;
712 
713  // cuts for each subdet_ector mimiking "Scheme B"
714  // double cutHB = 0.9, cutHE = 1.4, cutHO = 1.1, cutHFL = 1.2, cutHFS = 1.8;
715 
716  // energy in HCAL
717  double eHcal = 0.;
718  double eHcalCone = 0.;
719  double eHcalConeHB = 0.;
720  double eHcalConeHE = 0.;
721  double eHcalConeHO = 0.;
722  double eHcalConeHF = 0.;
723  double eHcalConeHFL = 0.;
724  double eHcalConeHFS = 0.;
725  // Total numbet of RecHits in HCAL, in the cone, above 1 GeV theshold
726  int nrechits = 0;
727  int nrechitsCone = 0;
728  int nrechitsThresh = 0;
729 
730  // energy in ECAL
731  double eEcal = 0.;
732  double eEcalB = 0.;
733  double eEcalE = 0.;
734  double eEcalCone = 0.;
735  int numrechitsEcal = 0;
736 
737  // MC info
738  double phi_MC = -999999.; // phi of initial particle from HepMC
739  double eta_MC = -999999.; // eta of initial particle from HepMC
740 
741  // HCAL energy around MC eta-phi at all depths;
742  double partR = 0.3;
743  double ehcal_coneMC_1 = 0.;
744  double ehcal_coneMC_2 = 0.;
745  double ehcal_coneMC_3 = 0.;
746  double ehcal_coneMC_4 = 0.;
747 
748  // Cone size for serach of the hottest HCAL cell around MC
749  double searchR = 1.0;
750  double eps = 0.001;
751 
752  // Single particle samples: actual eta-phi position of cluster around
753  // hottest cell
754  double etaHot = 99999.;
755  double phiHot = 99999.;
756 
757  // MC information
758 
759  // std::cout << "*** 1" << std::endl;
760 
761 
762  if(imc != 0) {
763 
765  ev.getByToken(tok_evt_,evtMC); // generator in late 310_preX
766  if (!evtMC.isValid()) {
767  std::cout << "no HepMCProduct found" << std::endl;
768  } else {
769  // std::cout << "*** source HepMCProduct found"<< std::endl;
770  }
771 
772  // MC particle with highest pt is taken as a direction reference
773  double maxPt = -99999.;
774  int npart = 0;
775  const HepMC::GenEvent * myGenEvent = evtMC->GetEvent();
776  for ( HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin();
777  p != myGenEvent->particles_end(); ++p ) {
778  double phip = (*p)->momentum().phi();
779  double etap = (*p)->momentum().eta();
780  // phi_MC = phip;
781  // eta_MC = etap;
782  double pt = (*p)->momentum().perp();
783  if(pt > maxPt) { npart++; maxPt = pt; phi_MC = phip; eta_MC = etap; }
784  }
785  // std::cout << "*** Max pT = " << maxPt << std::endl;
786 
787  }
788 
789  // std::cout << "*** 2" << std::endl;
790  // previously was: c.get<IdealGeometryRecord>().get (geometry);
791  c.get<CaloGeometryRecord>().get (geometry);
792 
793  // HCAL channel status map ****************************************
795  c.get<HcalChannelQualityRcd>().get( "withTopo", hcalChStatus );
796 
797  theHcalChStatus = hcalChStatus.product();
798 
799  // Assignment of severity levels **********************************
800  edm::ESHandle<HcalSeverityLevelComputer> hcalSevLvlComputerHndl;
801  c.get<HcalSeverityLevelComputerRcd>().get(hcalSevLvlComputerHndl);
802  theHcalSevLvlComputer = hcalSevLvlComputerHndl.product();
803 
804  // Fill working vectors of HCAL RecHits quantities (all of these are drawn)
805  fillRecHitsTmp(subdet_, ev);
806 
807  // HB
808  if( subdet_ ==5 || subdet_ == 1 ){
809  for(unsigned int iv=0; iv<hcalHBSevLvlVec.size(); iv++){
811  }
812  }
813  // HE
814  if( subdet_ ==5 || subdet_ == 2 ){
815  for(unsigned int iv=0; iv<hcalHESevLvlVec.size(); iv++){
817  }
818  }
819  // HO
820  if( subdet_ ==5 || subdet_ == 3 ){
821  for(unsigned int iv=0; iv<hcalHOSevLvlVec.size(); iv++){
823  }
824  }
825  // HF
826  if( subdet_ ==5 || subdet_ == 4 ){
827  for(unsigned int iv=0; iv<hcalHFSevLvlVec.size(); iv++){
829  }
830  }
831 
832  // std::cout << "*** 3" << std::endl;
833 
834 
835  //===========================================================================
836  // IN ALL other CASES : ieta-iphi maps
837  //===========================================================================
838 
839  // ECAL
840  if(ecalselector_ == "yes" && (subdet_ == 1 || subdet_ == 2 || subdet_ == 5)) {
842 
843 
844  ev.getByToken(tok_EB_, rhitEB);
845 
846  EcalRecHitCollection::const_iterator RecHit = rhitEB.product()->begin();
847  EcalRecHitCollection::const_iterator RecHitEnd = rhitEB.product()->end();
848 
849  for (; RecHit != RecHitEnd ; ++RecHit) {
850  EBDetId EBid = EBDetId(RecHit->id());
851 
852  const CaloCellGeometry* cellGeometry =
853  geometry->getSubdetectorGeometry (EBid)->getGeometry (EBid) ;
854  double eta = cellGeometry->getPosition ().eta () ;
855  double phi = cellGeometry->getPosition ().phi () ;
856  double en = RecHit->energy();
857  eEcal += en;
858  eEcalB += en;
859 
860  if (useAllHistos_) map_ecal->Fill(eta, phi, en);
861 
862  double r = dR(eta_MC, phi_MC, eta, phi);
863  if( r < partR) {
864  eEcalCone += en;
865  numrechitsEcal++;
866  }
867  }
868 
869 
871 
872  ev.getByToken(tok_EE_, rhitEE);
873 
874  RecHit = rhitEE.product()->begin();
875  RecHitEnd = rhitEE.product()->end();
876 
877  for (; RecHit != RecHitEnd ; ++RecHit) {
878  EEDetId EEid = EEDetId(RecHit->id());
879 
880  const CaloCellGeometry* cellGeometry =
881  geometry->getSubdetectorGeometry (EEid)->getGeometry (EEid) ;
882  double eta = cellGeometry->getPosition ().eta () ;
883  double phi = cellGeometry->getPosition ().phi () ;
884  double en = RecHit->energy();
885  eEcal += en;
886  eEcalE += en;
887 
888  if (useAllHistos_) map_ecal->Fill(eta, phi, en);
889 
890  double r = dR(eta_MC, phi_MC, eta, phi);
891  if( r < partR) {
892  eEcalCone += en;
893  numrechitsEcal++;
894  }
895  }
896  } // end of ECAL selection
897 
898 
899  // std::cout << "*** 4" << std::endl;
900 
901 
902  // Counting, including ZS items
903  // Filling HCAL maps ----------------------------------------------------
904  double maxE = -99999.;
905 
906  int nhb1 = 0;
907  int nhb2 = 0;
908  int nhe1 = 0;
909  int nhe2 = 0;
910  int nhe3 = 0;
911  int nho = 0;
912  int nhf1 = 0;
913  int nhf2 = 0;
914 
915  for (unsigned int i = 0; i < cen.size(); i++) {
916 
917  int sub = csub[i];
918  int depth = cdepth[i];
919  int ieta = cieta[i];
920  int iphi = ciphi[i];
921  double en = cen[i];
922  double eta = ceta[i];
923  double phi = cphi[i];
924  uint32_t stwd = cstwd[i];
925  uint32_t auxstwd = cauxstwd[i];
926  // double z = cz[i];
927 
928  int index = ieta * 72 + iphi; // for sequential histos
929 
930  /*
931  std::cout << "*** point 4-1" << " ieta, iphi, depth, sub = "
932  << ieta << ", " << iphi << ", " << depth << ", " << sub
933  << std::endl;
934  */
935 
936 
937  if( sub == 1 && depth == 1) nhb1++;
938  if( sub == 1 && depth == 2) nhb2++;
939  if( sub == 2 && depth == 1) nhe1++;
940  if( sub == 2 && depth == 2) nhe2++;
941  if( sub == 2 && depth == 3) nhe3++;
942  if( sub == 3 && depth == 4) nho++;
943  if( sub == 4 && depth == 1) nhf1++;
944  if( sub == 4 && depth == 2) nhf2++;
945 
946  if( subdet_ == 6) { // ZS specific
947  if( en < emap_min[ieta+41][iphi][depth-1][sub-1] )
948  emap_min[ieta+41][iphi][depth-1][sub-1] = en;
949  }
950 
951  double emin = 1.;
952  if(fabs(eta) > 3.) emin = 5.;
953 
954  double r = dR(eta_MC, phi_MC, eta, phi);
955  if( r < searchR ) { // search for hottest cell in a big cone
956  if(maxE < en && en > emin) {
957  maxE = en;
958  etaHot = eta;
959  phiHot = phi;
960  }
961  }
962 
963  /*
964  if(ieta == 27 ) {
965  std::cout << "*** ieta=28, iphi = " << iphi << " det = "
966  << sub << " depth = " << depth << std::endl;
967  }
968  */
969 
970  if( subdet_ != 6) {
971 
972  // std::cout << "*** 4-1" << std::endl;
973  //The emean_vs_ieta histos are drawn as well as the e_maps
974 
975 
976  // to distinguish HE and HF
977  if( depth == 1 || depth == 2 ) {
978  int ieta1 = ieta;
979  if(sub == 4) {
980  if (ieta1 < 0) ieta1--;
981  else ieta1++;
982  }
983  if (depth == 1) emap_depth1->Fill(double(ieta1), double(iphi), en);
984  if (depth == 2) emap_depth2->Fill(double(ieta1), double(iphi), en);
985  }
986 
987  if( depth == 3) emap_depth3->Fill(double(ieta), double(iphi), en);
988  if( depth == 4) emap_depth4->Fill(double(ieta), double(iphi), en);
989 
990  if (depth == 1 && sub == 1 ) {
991  emean_vs_ieta_HB1->Fill(double(ieta), en);
992  occupancy_map_HB1->Fill(double(ieta), double(iphi));
993  if(useAllHistos_){
994  emean_seqHB1->Fill(double(index), en);
995  }
996  }
997  if (depth == 2 && sub == 1) {
998  emean_vs_ieta_HB2->Fill(double(ieta), en);
999  occupancy_map_HB2->Fill(double(ieta), double(iphi));
1000  if(useAllHistos_){
1001  emean_seqHB2->Fill(double(index), en);
1002  }
1003  }
1004  if (depth == 1 && sub == 2) {
1005  emean_vs_ieta_HE1->Fill(double(ieta), en);
1006  occupancy_map_HE1->Fill(double(ieta), double(iphi));
1007  if(useAllHistos_){
1008  emean_seqHE1->Fill(double(index), en);
1009  }
1010  }
1011  if (depth == 2 && sub == 2) {
1012  emean_vs_ieta_HE2->Fill(double(ieta), en);
1013  occupancy_map_HE2->Fill(double(ieta), double(iphi));
1014  if(useAllHistos_){
1015  emean_seqHE2->Fill(double(index), en);
1016  }
1017  }
1018  if (depth == 3 && sub == 2) {
1019  emean_vs_ieta_HE3->Fill(double(ieta), en);
1020  occupancy_map_HE3->Fill(double(ieta), double(iphi));
1021  if(useAllHistos_){
1022  emean_seqHE3->Fill(double(index), en);
1023  }
1024  }
1025  if (depth == 4 ) {
1026  emean_vs_ieta_HO->Fill(double(ieta), en);
1027  occupancy_map_HO->Fill(double(ieta), double(iphi));
1028  if(useAllHistos_){
1029  emean_seqHO->Fill(double(index), en);
1030  }
1031  }
1032  if (depth == 1 && sub == 4) {
1033  emean_vs_ieta_HF1->Fill(double(ieta), en);
1034  occupancy_map_HF1->Fill(double(ieta), double(iphi));
1035  if(useAllHistos_){
1036  emean_seqHF1->Fill(double(index), en);
1037  }
1038  }
1039  if (depth == 2 && sub == 4) {
1040  emean_vs_ieta_HF2->Fill(double(ieta), en);
1041  occupancy_map_HF2->Fill(double(ieta), double(iphi));
1042  if(useAllHistos_){
1043  emean_seqHF2->Fill(double(index), en);
1044  }
1045  }
1046  }
1047 
1048 
1049  if( r < partR ) {
1050  if (depth == 1) ehcal_coneMC_1 += en;
1051  if (depth == 2) ehcal_coneMC_2 += en;
1052  if (depth == 3) ehcal_coneMC_3 += en;
1053  if (depth == 4) ehcal_coneMC_4 += en;
1054  }
1055 
1056  //32-bit status word
1057  uint32_t statadd;
1058  unsigned int isw67 = 0;
1059 
1060  //Status word correlation
1061  unsigned int sw27 = 27;
1062  unsigned int sw13 = 13;
1063 
1064  uint32_t statadd27 = 0x1<<sw27;
1065  uint32_t statadd13 = 0x1<<sw13;
1066 
1067  float status27 = 0;
1068  float status13 = 0;
1069 
1070  if(stwd & statadd27) status27 = 1;
1071  if(stwd & statadd13) status13 = 1;
1072 
1073  if (sub == 1){
1074  RecHit_StatusWordCorr_HB->Fill(status13, status27);
1075  } else if (sub == 2){
1076  RecHit_StatusWordCorr_HE->Fill(status13, status27);
1077  }
1078 
1079 
1080 
1081  for (unsigned int isw = 0; isw < 32; isw++){
1082  statadd = 0x1<<(isw);
1083  if (stwd & statadd){
1084  if (sub == 1) RecHit_StatusWord_HB->Fill(isw);
1085  else if (sub == 2) RecHit_StatusWord_HE->Fill(isw);
1086  else if (sub == 3) RecHit_StatusWord_HO->Fill(isw);
1087  else if (sub == 4){
1088  RecHit_StatusWord_HF->Fill(isw);
1089  if (isw == 6) isw67 += 1;
1090  if (isw == 7) isw67 += 2;
1091  }
1092  }
1093  }
1094  if (isw67 != 0 && useAllHistos_) RecHit_StatusWord_HF67->Fill(isw67); //This one is not drawn
1095 
1096  for (unsigned int isw =0; isw < 32; isw++){
1097  statadd = 0x1<<(isw);
1098  if( auxstwd & statadd ){
1099  if (sub == 1) RecHit_Aux_StatusWord_HB->Fill(isw);
1100  else if (sub == 2) RecHit_Aux_StatusWord_HE->Fill(isw);
1101  else if (sub == 3) RecHit_Aux_StatusWord_HO->Fill(isw);
1102  else if (sub == 4) RecHit_Aux_StatusWord_HF->Fill(isw);
1103  }
1104 
1105  }
1106 
1107  }
1108 
1109  // std::cout << "*** 4-2" << std::endl;
1110 
1111  if( subdet_ == 6 && useAllHistos_) { // ZS plots; not drawn
1112  ZS_nHB1->Fill(double(nhb1));
1113  ZS_nHB2->Fill(double(nhb2));
1114  ZS_nHE1->Fill(double(nhe1));
1115  ZS_nHE2->Fill(double(nhe2));
1116  ZS_nHE3->Fill(double(nhe3));
1117  ZS_nHO ->Fill(double(nho));
1118  ZS_nHF1->Fill(double(nhf1));
1119  ZS_nHF2->Fill(double(nhf2));
1120  }
1121  else{
1122  Nhb->Fill(double(nhb1 + nhb2));
1123  Nhe->Fill(double(nhe1 + nhe2 + nhe3));
1124  Nho->Fill(double(nho));
1125  Nhf->Fill(double(nhf1 + nhf2));
1126 
1127  //These are not drawn
1128  if(imc != 0 && useAllHistos_) {
1129  map_econe_depth1->Fill(eta_MC, phi_MC, ehcal_coneMC_1);
1130  map_econe_depth2->Fill(eta_MC, phi_MC, ehcal_coneMC_2);
1131  map_econe_depth3->Fill(eta_MC, phi_MC, ehcal_coneMC_3);
1132  map_econe_depth4->Fill(eta_MC, phi_MC, ehcal_coneMC_4);
1133  }
1134  }
1135 
1136  // std::cout << "*** 5" << std::endl;
1137 
1138 
1139  // NOISE =================================================================
1140  //Not drawn
1141  if (hcalselector_ == "noise" && useAllHistos_) {
1142  for (unsigned int i = 0; i < cen.size(); i++) {
1143 
1144  int sub = csub[i];
1145  int depth = cdepth[i];
1146  double en = cen[i];
1147 
1148  if (sub == 1) e_hb->Fill(en);
1149  if (sub == 2) e_he->Fill(en);
1150  if (sub == 3) e_ho->Fill(en);
1151  if (sub == 4) {
1152  if(depth == 1)
1153  e_hfl->Fill(en);
1154  else
1155  e_hfs->Fill(en);
1156  }
1157  }
1158  }
1159 
1160  //===========================================================================
1161  // SUBSYSTEMS,
1162  //===========================================================================
1163 
1164  else if ((subdet_ != 6) && (subdet_ != 0)) {
1165 
1166  // std::cout << "*** 6" << std::endl;
1167 
1168 
1169  double clusEta = 999.;
1170  double clusPhi = 999.;
1171  double clusEn = 0.;
1172 
1173  double HcalCone_d1 = 0.;
1174  double HcalCone_d2 = 0.;
1175  double HcalCone_d3 = 0.;
1176  double HcalCone_d4 = 0.;
1177  double HcalCone = 0.;
1178 
1179  int ietaMax1 = 9999;
1180  int ietaMax2 = 9999;
1181  int ietaMax3 = 9999;
1182  int ietaMax4 = 9999;
1183  int ietaMax = 9999;
1184  double enMax1 = -9999.;
1185  double enMax2 = -9999.;
1186  double enMax3 = -9999.;
1187  double enMax4 = -9999.;
1188  // double enMax = -9999.;
1189  double etaMax = 9999.;
1190 
1191  /*
1192  std::cout << "*** point 5-1" << " eta_MC, phi_MC etaHot, phiHot = "
1193  << eta_MC << ", " << phi_MC << " "
1194  << etaHot << ", " << phiHot
1195  << std::endl;
1196  */
1197 
1198  // CYCLE over cells ====================================================
1199 
1200  for (unsigned int i = 0; i < cen.size(); i++) {
1201  int sub = csub[i];
1202  int depth = cdepth[i];
1203  double eta = ceta[i];
1204  double phi = cphi[i];
1205  double en = cen[i];
1206  double t = ctime[i];
1207  int ieta = cieta[i];
1208 
1209  double rhot = dR(etaHot, phiHot, eta, phi);
1210  if(rhot < partR && en > 1.) {
1211  clusEta = (clusEta * clusEn + eta * en)/(clusEn + en);
1212  clusPhi = phi12(clusPhi, clusEn, phi, en);
1213  clusEn += en;
1214  }
1215 
1216  nrechits++;
1217  eHcal += en;
1218  if(en > 1. ) nrechitsThresh++;
1219 
1220  double r = dR(eta_MC, phi_MC, eta, phi);
1221  if( r < partR ){
1222  if(sub == 1) eHcalConeHB += en;
1223  if(sub == 2) eHcalConeHE += en;
1224  if(sub == 3) eHcalConeHO += en;
1225  if(sub == 4) {
1226  eHcalConeHF += en;
1227  if (depth == 1) eHcalConeHFL += en;
1228  else eHcalConeHFS += en;
1229  }
1230  eHcalCone += en;
1231  nrechitsCone++;
1232 
1233  // search for most energetic cell at the given depth in the cone
1234  if(depth == 1) {
1235  HcalCone_d1 += en;
1236  if(enMax1 < en) {
1237  enMax1 = en;
1238  ietaMax1 = ieta;
1239  }
1240  }
1241  if(depth == 2) {
1242  HcalCone_d2 += en;
1243  if(enMax2 < en) {
1244  enMax2 = en;
1245  ietaMax2 = ieta;
1246  }
1247  }
1248  if(depth == 3) {
1249  HcalCone_d3 += en;
1250  if(enMax3 < en) {
1251  enMax3 = en;
1252  ietaMax3 = ieta;
1253  }
1254  }
1255  if(depth == 4) {
1256  HcalCone_d4 += en;
1257  if(enMax4 < en) {
1258  enMax4 = en;
1259  ietaMax4 = ieta;
1260  }
1261  }
1262 
1263  if(depth != 4) {
1264  HcalCone += en;
1265  }
1266 
1267 
1268  // regardless of the depths (but excluding HO), just hottest cell
1269  /*
1270  if(depth != 4) {
1271  if(enMax < en) {
1272  enMax = en;
1273  ietaMax = ieta;
1274  }
1275  }
1276  */
1277 
1278  // alternative: ietamax -> closest to MC eta !!!
1279  float eta_diff = fabs(eta_MC - eta);
1280  if(eta_diff < etaMax) {
1281  etaMax = eta_diff;
1282  ietaMax = ieta;
1283  }
1284  }
1285 
1286  //The energy and overall timing histos are drawn while
1287  //the ones split by depth are not
1288  if(sub == 1 && (subdet_ == 1 || subdet_ == 5)) {
1289  meTimeHB->Fill(t);
1290  meRecHitsEnergyHB->Fill(en);
1291 
1292  meTE_Low_HB->Fill( en, t);
1293  meTE_HB->Fill( en, t);
1294  meTE_High_HB->Fill( en, t);
1295  meTEprofileHB_Low->Fill(en, t);
1296  meTEprofileHB->Fill(en, t);
1297  meTEprofileHB_High->Fill(en, t);
1298 
1299  if (useAllHistos_){
1300  if (depth == 1) meTE_HB1->Fill( en, t);
1301  else if (depth == 2) meTE_HB2->Fill( en, t);
1302  }
1303  }
1304  if(sub == 2 && (subdet_ == 2 || subdet_ == 5)) {
1305  meTimeHE->Fill(t);
1306  meRecHitsEnergyHE->Fill(en);
1307 
1308  meTE_Low_HE->Fill( en, t);
1309  meTE_HE->Fill( en, t);
1310  meTEprofileHE_Low->Fill(en, t);
1311  meTEprofileHE->Fill(en, t);
1312 
1313  if (useAllHistos_){
1314  if (depth == 1) meTE_HE1->Fill( en, t);
1315  else if (depth == 2) meTE_HE2->Fill( en, t);
1316  }
1317  }
1318  if(sub == 4 && (subdet_ == 4 || subdet_ == 5)) {
1319  meTimeHF->Fill(t);
1320  meRecHitsEnergyHF->Fill(en);
1321 
1322  meTE_Low_HF->Fill(en, t);
1323  meTE_HF->Fill(en, t);
1324  meTEprofileHF_Low->Fill(en, t);
1325  meTEprofileHF->Fill(en, t);
1326 
1327  if (useAllHistos_){
1328  if (depth == 1) meTE_HFL->Fill( en, t);
1329  else meTE_HFS->Fill( en, t);
1330  }
1331  }
1332  if(sub == 3 && (subdet_ == 3 || subdet_ == 5)) {
1333  meTimeHO->Fill(t);
1334  meRecHitsEnergyHO->Fill(en);
1335 
1336  meTE_HO->Fill( en, t);
1337  meTE_High_HO->Fill( en, t);
1338  meTEprofileHO->Fill(en, t);
1339  meTEprofileHO_High->Fill(en, t);
1340  }
1341  }
1342 
1343  if(imc != 0) {
1344  //Cone by depth are not drawn, the others are used for pion scan
1345  if (useAllHistos_){
1346  meEnConeEtaProfile_depth1->Fill(double(ietaMax1), HcalCone_d1);
1347  meEnConeEtaProfile_depth2->Fill(double(ietaMax2), HcalCone_d2);
1348  meEnConeEtaProfile_depth3->Fill(double(ietaMax3), HcalCone_d3);
1349  meEnConeEtaProfile_depth4->Fill(double(ietaMax4), HcalCone_d4);
1350  }
1351  meEnConeEtaProfile ->Fill(double(ietaMax), HcalCone); //
1352  meEnConeEtaProfile_E ->Fill(double(ietaMax), eEcalCone);
1353  meEnConeEtaProfile_EH ->Fill(double(ietaMax), HcalCone+eEcalCone);
1354  }
1355 
1356  // std::cout << "*** 7" << std::endl;
1357 
1358 
1359  // Single particle samples ONLY ! ======================================
1360  // Fill up some histos for "integrated" subsustems.
1361  // These are not drawn
1362  if(etype_ == 1 && useAllHistos_) {
1363 
1364  /*
1365  std::cout << "*** point 7-1" << " eta_MC, phi_MC clusEta, clusPhi = "
1366  << eta_MC << ", " << phi_MC << " "
1367  << clusEta << ", " << clusPhi
1368  << std::endl;
1369  */
1370 
1371  double phidev = dPhiWsign(clusPhi, phi_MC);
1372  meDeltaPhi->Fill(eta_MC, phidev);
1373  double etadev = clusEta - eta_MC;
1374  meDeltaEta->Fill(eta_MC, etadev);
1375 
1376  if(subdet_ == 1) {
1377  meSumRecHitsEnergyHB->Fill(eHcal);
1378  if(imc != 0) meSumRecHitsEnergyConeHB->Fill(eHcalConeHB);
1379  if(imc != 0) meNumRecHitsConeHB->Fill(double(nrechitsCone));
1380  meNumRecHitsThreshHB->Fill(double(nrechitsThresh));
1381  }
1382 
1383  if(subdet_ == 2) {
1384  meSumRecHitsEnergyHE->Fill(eHcal);
1385  if(imc != 0) meSumRecHitsEnergyConeHE->Fill(eHcalConeHE);
1386  if(imc != 0) meNumRecHitsConeHE->Fill(double(nrechitsCone));
1387  meNumRecHitsThreshHE->Fill(double(nrechitsThresh));
1388  }
1389 
1390  if(subdet_ == 3) {
1391  meSumRecHitsEnergyHO->Fill(eHcal);
1392  if(imc != 0) meSumRecHitsEnergyConeHO->Fill(eHcalConeHO);
1393  if(imc != 0) meNumRecHitsConeHO->Fill(double(nrechitsCone));
1394  meNumRecHitsThreshHO->Fill(double(nrechitsThresh));
1395  }
1396 
1397  if(subdet_ == 4) {
1398  if(eHcalConeHF > eps ) {
1399  meSumRecHitsEnergyHF ->Fill(eHcal);
1400  if(imc != 0) {
1401  meSumRecHitsEnergyConeHF ->Fill(eHcalConeHF);
1402  meNumRecHitsConeHF->Fill(double(nrechitsCone));
1403  meSumRecHitsEnergyConeHFL ->Fill(eHcalConeHFL);
1404  meSumRecHitsEnergyConeHFS ->Fill(eHcalConeHFS);
1405  }
1406  }
1407  }
1408 
1409  // std::cout << "*** 8" << std::endl;
1410 
1411 
1412  // Also combine with ECAL if needed
1413  if(subdet_ == 1 && ecalselector_ == "yes") {
1414 
1415  /*
1416  std::cout << "*** point 8-1"
1417  << " eEcalB " << eEcalB << " eHcal " << eHcal
1418  << " eEcalCone " << eEcalCone << " eHcalCone "
1419  << eHcalCone
1420  << " numrechitsEcal " << numrechitsEcal
1421  << std::endl;
1422 
1423  */
1424 
1425  meEcalHcalEnergyHB->Fill(eEcalB+eHcal);
1426  meEcalHcalEnergyConeHB->Fill(eEcalCone+eHcalCone);
1427  meNumEcalRecHitsConeHB->Fill(double(numrechitsEcal));
1428 
1429  }
1430 
1431  if(subdet_ == 2 && ecalselector_ == "yes"){
1432 
1433  /*
1434  std::cout << "*** point 8-2a"
1435  << " eEcalE " << eEcalE << " eHcal " << eHcal
1436  << " eEcalCone " << eEcalCone << " eHcalCone "
1437  << eHcalCone
1438  << " numrechitsEcal " << numrechitsEcal
1439  << std::endl;
1440  */
1441 
1442  meEcalHcalEnergyHE->Fill(eEcalE+eHcal);
1443  if(imc != 0) meEcalHcalEnergyConeHE->Fill(eEcalCone+eHcalCone);
1444  if(imc != 0) meNumEcalRecHitsConeHE->Fill(double(numrechitsEcal));
1445  }
1446 
1447  // Banana plots finally
1448  if(imc != 0) {
1449  if(subdet_ == 1 && ecalselector_ == "yes")
1450  meEnergyHcalVsEcalHB -> Fill(eEcalCone,eHcalCone);
1451  if(subdet_ == 2 && ecalselector_ == "yes")
1452  meEnergyHcalVsEcalHE -> Fill(eEcalCone,eHcalCone);
1453  }
1454  }
1455  }
1456  // std::cout << "*** 9" << std::endl;
1457 
1458 
1459  //===========================================================================
1460  // Getting SimHits
1461  //===========================================================================
1462 
1463  if(subdet_ > 0 && subdet_ < 6 && imc !=0 && !famos_ ) { // not noise
1464 
1465  double maxES = -9999.;
1466  double etaHotS = 1000.;
1467  double phiHotS = 1000.;
1468 
1470  ev.getByToken(tok_hh_,hcalHits);
1471  const PCaloHitContainer * SimHitResult = hcalHits.product () ;
1472 
1473  double enSimHits = 0.;
1474  double enSimHitsHB = 0.;
1475  double enSimHitsHE = 0.;
1476  double enSimHitsHO = 0.;
1477  double enSimHitsHF = 0.;
1478  double enSimHitsHFL = 0.;
1479  double enSimHitsHFS = 0.;
1480  // sum of SimHits in the cone
1481 
1482  for (std::vector<PCaloHit>::const_iterator SimHits = SimHitResult->begin () ; SimHits != SimHitResult->end(); ++SimHits) {
1483  HcalDetId cell(SimHits->id());
1484  int sub = cell.subdet();
1485  const CaloCellGeometry* cellGeometry =
1486  geometry->getSubdetectorGeometry (cell)->getGeometry (cell) ;
1487  double etaS = cellGeometry->getPosition().eta () ;
1488  double phiS = cellGeometry->getPosition().phi () ;
1489  double en = SimHits->energy();
1490 
1491  double emin = 0.01;
1492  if(fabs(etaS) > 3.) emin = 1.;
1493 
1494  double r = dR(eta_MC, phi_MC, etaS, phiS);
1495  if( r < searchR ) { // search for hottest cell in a big cone
1496  if(maxES < en && en > emin ) {
1497  maxES = en;
1498  etaHotS = etaS;
1499  phiHotS = phiS;
1500  }
1501  }
1502 
1503  if ( r < partR ){ // just energy in the small cone
1504  enSimHits += en;
1505  if(sub == 1) enSimHitsHB += en;
1506  if(sub == 2) enSimHitsHE += en;
1507  if(sub == 3) enSimHitsHO += en;
1508  if(sub == 4) {
1509  enSimHitsHF += en;
1510  int depth = cell.depth();
1511  if(depth == 1) enSimHitsHFL += en;
1512  else enSimHitsHFS += en;
1513  }
1514  }
1515  }
1516 
1517 
1518  // Second look over SimHits: cluster finding
1519 
1520  double clusEta = 999.;
1521  double clusPhi = 999.;
1522  double clusEn = 0.;
1523 
1524  for (std::vector<PCaloHit>::const_iterator SimHits = SimHitResult->begin () ; SimHits != SimHitResult->end(); ++SimHits) {
1525  HcalDetId cell(SimHits->id());
1526 
1527  const CaloCellGeometry* cellGeometry =
1528  geometry->getSubdetectorGeometry (cell)->getGeometry (cell) ;
1529  double etaS = cellGeometry->getPosition().eta () ;
1530  double phiS = cellGeometry->getPosition().phi () ;
1531  double en = SimHits->energy();
1532 
1533  double emin = 0.01;
1534  if(fabs(etaS) > 3.) emin = 1.;
1535 
1536  double rhot = dR(etaHotS, phiHotS, etaS, phiS);
1537  if(rhot < partR && en > emin) {
1538  clusEta = (clusEta * clusEn + etaS * en)/(clusEn + en);
1539  clusPhi = phi12(clusPhi, clusEn, phiS, en);
1540  clusEn += en;
1541  }
1542  }
1543 
1544  // SimHits cluster deviation from MC (eta, phi)
1545  // These are not drawn
1546  if (useAllHistos_){
1547  if(etype_ == 1) {
1548  double phidev = dPhiWsign(clusPhi, phi_MC);
1549  meDeltaPhiS->Fill(eta_MC, phidev);
1550  double etadev = clusEta - eta_MC;
1551  meDeltaEtaS->Fill(eta_MC, etadev);
1552  }
1553  // Now some histos with SimHits
1554 
1555  if(subdet_ == 4 || subdet_ == 5) {
1556  if(eHcalConeHF > eps) {
1557  meRecHitSimHitHF->Fill( enSimHitsHF, eHcalConeHF );
1558  meRecHitSimHitProfileHF->Fill( enSimHitsHF, eHcalConeHF);
1559 
1560  meRecHitSimHitHFL->Fill( enSimHitsHFL, eHcalConeHFL );
1561  meRecHitSimHitProfileHFL->Fill( enSimHitsHFL, eHcalConeHFL);
1562  meRecHitSimHitHFS->Fill( enSimHitsHFS, eHcalConeHFS );
1563  meRecHitSimHitProfileHFS->Fill( enSimHitsHFS, eHcalConeHFS);
1564  }
1565  }
1566  if(subdet_ == 1 || subdet_ == 5) {
1567  meRecHitSimHitHB->Fill( enSimHitsHB,eHcalConeHB );
1568  meRecHitSimHitProfileHB->Fill( enSimHitsHB,eHcalConeHB);
1569  }
1570  if(subdet_ == 2 || subdet_ == 5) {
1571  meRecHitSimHitHE->Fill( enSimHitsHE,eHcalConeHE );
1572  meRecHitSimHitProfileHE->Fill( enSimHitsHE,eHcalConeHE);
1573  }
1574  if(subdet_ == 3 || subdet_ == 5) {
1575  meRecHitSimHitHO->Fill( enSimHitsHO,eHcalConeHO );
1576  meRecHitSimHitProfileHO->Fill( enSimHitsHO,eHcalConeHO);
1577  }
1578  }
1579  }
1580 
1581  nevtot++;
1582 }
1583 
1584 
1587 
1588  using namespace edm;
1589 
1590 
1591  // initialize data vectors
1592  csub.clear();
1593  cen.clear();
1594  ceta.clear();
1595  cphi.clear();
1596  ctime.clear();
1597  cieta.clear();
1598  ciphi.clear();
1599  cdepth.clear();
1600  cz.clear();
1601  cstwd.clear();
1602  cauxstwd.clear();
1603  hcalHBSevLvlVec.clear();
1604  hcalHESevLvlVec.clear();
1605  hcalHFSevLvlVec.clear();
1606  hcalHOSevLvlVec.clear();
1607 
1608  if( subdet_ == 1 || subdet_ == 2 || subdet_ == 5 || subdet_ == 6 || subdet_ == 0) {
1609 
1610  //HBHE
1612  ev.getByToken(tok_hbhe_, hbhecoll);
1613 
1614  for (HBHERecHitCollection::const_iterator j=hbhecoll->begin(); j != hbhecoll->end(); j++) {
1615  HcalDetId cell(j->id());
1616  const CaloCellGeometry* cellGeometry =
1617  geometry->getSubdetectorGeometry (cell)->getGeometry (cell) ;
1618  double eta = cellGeometry->getPosition().eta () ;
1619  double phi = cellGeometry->getPosition().phi () ;
1620  double zc = cellGeometry->getPosition().z ();
1621  int sub = cell.subdet();
1622  int depth = cell.depth();
1623  int inteta = cell.ieta();
1624  if(inteta > 0) inteta -= 1;
1625  int intphi = cell.iphi()-1;
1626  double en = j->energy();
1627  double t = j->time();
1628  int stwd = j->flags();
1629  int auxstwd = j->aux();
1630 
1631  int severityLevel = hcalSevLvl( (CaloRecHit*) &*j );
1632  if( cell.subdet()==HcalBarrel ){
1633  hcalHBSevLvlVec.push_back(severityLevel);
1634  }else if (cell.subdet()==HcalEndcap ){
1635  hcalHESevLvlVec.push_back(severityLevel);
1636  }
1637 
1638  if((iz > 0 && eta > 0.) || (iz < 0 && eta <0.) || iz == 0) {
1639 
1640  csub.push_back(sub);
1641  cen.push_back(en);
1642  ceta.push_back(eta);
1643  cphi.push_back(phi);
1644  ctime.push_back(t);
1645  cieta.push_back(inteta);
1646  ciphi.push_back(intphi);
1647  cdepth.push_back(depth);
1648  cz.push_back(zc);
1649  cstwd.push_back(stwd);
1650  cauxstwd.push_back(auxstwd);
1651  }
1652  }
1653 
1654  }
1655 
1656  if( subdet_ == 4 || subdet_ == 5 || subdet_ == 6 || subdet_ == 0) {
1657 
1658  //HF
1660  ev.getByToken(tok_hf_, hfcoll);
1661 
1662  for (HFRecHitCollection::const_iterator j = hfcoll->begin(); j != hfcoll->end(); j++) {
1663  HcalDetId cell(j->id());
1664  const CaloCellGeometry* cellGeometry =
1665  geometry->getSubdetectorGeometry (cell)->getGeometry (cell) ;
1666  double eta = cellGeometry->getPosition().eta () ;
1667  double phi = cellGeometry->getPosition().phi () ;
1668  double zc = cellGeometry->getPosition().z ();
1669  int sub = cell.subdet();
1670  int depth = cell.depth();
1671  int inteta = cell.ieta();
1672  if(inteta > 0) inteta -= 1;
1673  int intphi = cell.iphi()-1;
1674  double en = j->energy();
1675  double t = j->time();
1676  int stwd = j->flags();
1677  int auxstwd = j->aux();
1678 
1679  int severityLevel = hcalSevLvl( (CaloRecHit*) &*j );
1680  if( cell.subdet()==HcalForward ){
1681  hcalHFSevLvlVec.push_back(severityLevel);
1682  }
1683 
1684  if((iz > 0 && eta > 0.) || (iz < 0 && eta <0.) || iz == 0) {
1685 
1686  csub.push_back(sub);
1687  cen.push_back(en);
1688  ceta.push_back(eta);
1689  cphi.push_back(phi);
1690  ctime.push_back(t);
1691  cieta.push_back(inteta);
1692  ciphi.push_back(intphi);
1693  cdepth.push_back(depth);
1694  cz.push_back(zc);
1695  cstwd.push_back(stwd);
1696  cauxstwd.push_back(auxstwd);
1697  }
1698  }
1699  }
1700 
1701  //HO
1702  if( subdet_ == 3 || subdet_ == 5 || subdet_ == 6 || subdet_ == 0) {
1703 
1705  ev.getByToken(tok_ho_, hocoll);
1706 
1707  for (HORecHitCollection::const_iterator j = hocoll->begin(); j != hocoll->end(); j++) {
1708  HcalDetId cell(j->id());
1709  const CaloCellGeometry* cellGeometry =
1710  geometry->getSubdetectorGeometry (cell)->getGeometry (cell) ;
1711  double eta = cellGeometry->getPosition().eta () ;
1712  double phi = cellGeometry->getPosition().phi () ;
1713  double zc = cellGeometry->getPosition().z ();
1714  int sub = cell.subdet();
1715  int depth = cell.depth();
1716  int inteta = cell.ieta();
1717  if(inteta > 0) inteta -= 1;
1718  int intphi = cell.iphi()-1;
1719  double t = j->time();
1720  double en = j->energy();
1721  int stwd = j->flags();
1722  int auxstwd = j->aux();
1723 
1724  int severityLevel = hcalSevLvl( (CaloRecHit*) &*j );
1725  if( cell.subdet()==HcalOuter ){
1726  hcalHOSevLvlVec.push_back(severityLevel);
1727  }
1728 
1729  if((iz > 0 && eta > 0.) || (iz < 0 && eta <0.) || iz == 0) {
1730  csub.push_back(sub);
1731  cen.push_back(en);
1732  ceta.push_back(eta);
1733  cphi.push_back(phi);
1734  ctime.push_back(t);
1735  cieta.push_back(inteta);
1736  ciphi.push_back(intphi);
1737  cdepth.push_back(depth);
1738  cz.push_back(zc);
1739  cstwd.push_back(stwd);
1740  cauxstwd.push_back(auxstwd);
1741  }
1742  }
1743  }
1744 }
1745 
1746 double HcalRecHitsValidation::dR(double eta1, double phi1, double eta2, double phi2) {
1747  double PI = 3.1415926535898;
1748  double deltaphi= phi1 - phi2;
1749  if( phi2 > phi1 ) { deltaphi= phi2 - phi1;}
1750  if(deltaphi > PI) { deltaphi = 2.*PI - deltaphi;}
1751  double deltaeta = eta2 - eta1;
1752  double tmp = sqrt(deltaeta* deltaeta + deltaphi*deltaphi);
1753  return tmp;
1754 }
1755 
1756 double HcalRecHitsValidation::phi12(double phi1, double en1, double phi2, double en2) {
1757  // weighted mean value of phi1 and phi2
1758 
1759  double tmp;
1760  double PI = 3.1415926535898;
1761  double a1 = phi1; double a2 = phi2;
1762 
1763  if( a1 > 0.5*PI && a2 < 0.) a2 += 2*PI;
1764  if( a2 > 0.5*PI && a1 < 0.) a1 += 2*PI;
1765  tmp = (a1 * en1 + a2 * en2)/(en1 + en2);
1766  if(tmp > PI) tmp -= 2.*PI;
1767 
1768  return tmp;
1769 
1770 }
1771 
1772 double HcalRecHitsValidation::dPhiWsign(double phi1, double phi2) {
1773  // clockwise phi2 w.r.t phi1 means "+" phi distance
1774  // anti-clockwise phi2 w.r.t phi1 means "-" phi distance
1775 
1776  double PI = 3.1415926535898;
1777  double a1 = phi1; double a2 = phi2;
1778  double tmp = a2 - a1;
1779  if( a1*a2 < 0.) {
1780  if(a1 > 0.5 * PI) tmp += 2.*PI;
1781  if(a2 > 0.5 * PI) tmp -= 2.*PI;
1782  }
1783  return tmp;
1784 
1785 }
1786 
1788 
1789  const DetId id = hit->detid();
1790 
1791  const uint32_t recHitFlag = hit->flags();
1792  const uint32_t dbStatusFlag = theHcalChStatus->getValues(id)->getValue();
1793 
1794  int severityLevel = theHcalSevLvlComputer->getSeverityLevel(id, recHitFlag, dbStatusFlag);
1795 
1796  return severityLevel;
1797 
1798 }
1799 
1801 
MonitorElement * emean_vs_ieta_HB2
MonitorElement * RMS_vs_ieta_HB1
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
MonitorElement * occupancy_vs_ieta_HO
const HcalSeverityLevelComputer * theHcalSevLvlComputer
MonitorElement * meSumRecHitsEnergyConeHE
MonitorElement * occupancy_seqHB1
const HcalChannelQuality * theHcalChStatus
std::vector< PCaloHit > PCaloHitContainer
MonitorElement * occupancy_map_HB1
MonitorElement * occupancy_vs_ieta_HE3
edm::EDGetTokenT< edm::PCaloHitContainer > tok_hh_
int ib
Definition: cuy.py:660
MonitorElement * meRecHitsEnergyHF
MonitorElement * occupancy_vs_ieta_HF1
MonitorElement * RMS_vs_ieta_HF1
HcalSubdetector subdet() const
get the subdetector
Definition: HcalDetId.h:45
#define PI
double phi12(double phi1, double en1, double phi2, double en2)
MonitorElement * bookProfile(Args &&...args)
Definition: DQMStore.h:157
const DetId & detid() const
Definition: CaloRecHit.h:20
std::vector< int > hcalHBSevLvlVec
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:464
MonitorElement * meNumRecHitsConeHE
MonitorElement * meSumRecHitsEnergyConeHB
MonitorElement * meTEprofileHE_Low
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
std::vector< uint32_t > cstwd
MonitorElement * RMS_vs_ieta_HE3
MonitorElement * RecHit_Aux_StatusWord_HB
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
MonitorElement * RMS_vs_ieta_HE1
MonitorElement * emean_vs_ieta_HB1
edm::EDGetTokenT< HFRecHitCollection > tok_hf_
std::vector< EcalRecHit >::const_iterator const_iterator
MonitorElement * occupancy_seqHO
MonitorElement * occupancy_seqHE3
MonitorElement * meEnConeEtaProfile
double npart
Definition: HydjetWrapper.h:49
bool ev
MonitorElement * occupancy_map_HE3
double dPhiWsign(double phi1, double phi2)
MonitorElement * meRecHitSimHitProfileHF
const Item * getValues(DetId fId, bool throwOnFail=true) const
MonitorElement * meEcalHcalEnergyHB
int hcalSevLvl(const CaloRecHit *hit)
MonitorElement * meRecHitSimHitProfileHFS
MonitorElement * meEnergyHcalVsEcalHE
MonitorElement * occupancy_map_HF1
MonitorElement * meEcalHcalEnergyConeHE
MonitorElement * meRecHitsEnergyHB
MonitorElement * meSumRecHitsEnergyHF
double emap_min[82][72][4][4]
MonitorElement * meNumRecHitsConeHO
std::vector< double > ceta
MonitorElement * occupancy_seqHF2
MonitorElement * meNumRecHitsConeHB
MonitorElement * RMS_vs_ieta_HB2
MonitorElement * meRecHitSimHitHF
MonitorElement * meRecHitsEnergyHE
MonitorElement * occupancy_map_HE1
MonitorElement * meRecHitSimHitProfileHE
MonitorElement * meSumRecHitsEnergyHB
edm::EDGetTokenT< HORecHitCollection > tok_ho_
void Fill(long long x)
MonitorElement * meEnConeEtaProfile_E
MonitorElement * meRecHitSimHitHB
MonitorElement * meRecHitSimHitHFS
MonitorElement * RMS_vs_ieta_HO
MonitorElement * meSumRecHitsEnergyConeHFS
MonitorElement * RecHit_Aux_StatusWord_HO
std::vector< double > cphi
MonitorElement * meNumEcalRecHitsConeHB
MonitorElement * meNumRecHitsThreshHB
MonitorElement * meSumRecHitsEnergyHO
std::vector< int > hcalHESevLvlVec
MonitorElement * meEnergyHcalVsEcalHB
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
MonitorElement * RecHit_StatusWord_HF
T sqrt(T t)
Definition: SSEVec.h:48
MonitorElement * RecHit_Aux_StatusWord_HE
edm::ESHandle< CaloGeometry > geometry
MonitorElement * meSumRecHitsEnergyConeHO
MonitorElement * meTEprofileHB
uint32_t flags() const
Definition: CaloRecHit.h:21
MonitorElement * map_econe_depth1
MonitorElement * book1D(Args &&...args)
Definition: DQMStore.h:115
virtual void fillRecHitsTmp(int subdet_, edm::Event const &ev)
int j
Definition: DBlmapReader.cc:9
MonitorElement * map_econe_depth4
double dR(double eta1, double phi1, double eta2, double phi2)
MonitorElement * emean_vs_ieta_HE3
MonitorElement * meNumRecHitsThreshHO
MonitorElement * meEcalHcalEnergyConeHB
MonitorElement * meNumRecHitsThreshHE
MonitorElement * RecHit_StatusWordCorr_HB
edm::EDGetTokenT< HBHERecHitCollection > tok_hbhe_
std::vector< int > hcalHFSevLvlVec
MonitorElement * occupancy_map_HE2
MonitorElement * occupancy_vs_ieta_HB2
bool isValid() const
Definition: HandleBase.h:75
MonitorElement * meSumRecHitsEnergyConeHF
MonitorElement * RecHit_Aux_StatusWord_HF
MonitorElement * meRecHitSimHitHO
MonitorElement * meRecHitSimHitProfileHB
MonitorElement * RecHit_StatusWord_HO
MonitorElement * meEnConeEtaProfile_depth1
tuple conf
Definition: dbtoconf.py:185
MonitorElement * emean_vs_ieta_HF2
virtual void analyze(edm::Event const &ev, edm::EventSetup const &c)
virtual void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &)
MonitorElement * occupancy_map_HO
MonitorElement * occupancy_vs_ieta_HE2
MonitorElement * meTEprofileHF_Low
edm::EDGetTokenT< EBRecHitCollection > tok_EB_
MonitorElement * occupancy_map_HF2
MonitorElement * RecHit_StatusWord_HF67
MonitorElement * meTEprofileHO_High
MonitorElement * meTEprofileHE
Definition: DetId.h:18
MonitorElement * occupancy_seqHB2
MonitorElement * occupancy_map_HB2
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:274
MonitorElement * RMS_vs_ieta_HE2
T const * product() const
Definition: Handle.h:81
MonitorElement * meNumRecHitsConeHF
MonitorElement * book2D(Args &&...args)
Definition: DQMStore.h:133
MonitorElement * meTEprofileHF
MonitorElement * meEnConeEtaProfile_depth3
MonitorElement * emean_vs_ieta_HE1
MonitorElement * meRecHitSimHitProfileHO
std::vector< int > hcalHOSevLvlVec
MonitorElement * emean_vs_ieta_HE2
MonitorElement * occupancy_vs_ieta_HF2
edm::EDGetTokenT< EERecHitCollection > tok_EE_
const T & get() const
Definition: EventSetup.h:56
std::vector< double > ctime
MonitorElement * map_econe_depth2
int getSeverityLevel(const DetId &myid, const uint32_t &myflag, const uint32_t &mystatus) const
MonitorElement * RecHit_StatusWord_HE
MonitorElement * meEnConeEtaProfile_EH
MonitorElement * meRecHitSimHitProfileHFL
std::vector< double > cz
MonitorElement * meEnConeEtaProfile_depth2
T eta() const
Definition: PV3DBase.h:76
MonitorElement * RMS_vs_ieta_HF2
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
MonitorElement * meRecHitsEnergyHO
std::vector< uint32_t > cauxstwd
MonitorElement * emean_vs_ieta_HO
MonitorElement * occupancy_seqHF1
MonitorElement * meSumRecHitsEnergyHE
MonitorElement * meRecHitSimHitHFL
MonitorElement * meSumRecHitsEnergyConeHFL
MonitorElement * meNumEcalRecHitsConeHE
MonitorElement * meEcalHcalEnergyHE
MonitorElement * emean_vs_ieta_HF1
MonitorElement * meTEprofileHO
tuple cout
Definition: gather_cfg.py:121
edm::EDGetTokenT< edm::HepMCProduct > tok_evt_
MonitorElement * meTEprofileHB_Low
MonitorElement * RecHit_StatusWordCorr_HE
MonitorElement * occupancy_vs_ieta_HB1
uint32_t getValue() const
const GlobalPoint & getPosition() const
Returns the position of reference for this cell.
MonitorElement * occupancy_seqHE2
MonitorElement * meRecHitSimHitHE
HcalRecHitsValidation(edm::ParameterSet const &conf)
MonitorElement * map_econe_depth3
MonitorElement * occupancy_seqHE1
MonitorElement * meTEprofileHB_High
std::vector< double > cen
Definition: Run.h:43
MonitorElement * meEnConeEtaProfile_depth4
MonitorElement * occupancy_vs_ieta_HE1
MonitorElement * RecHit_StatusWord_HB