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