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